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Brief Summary of Work Performed Since August 1 . 1993 : 

We have systematically investigated the microprocesses occurring when a magnetic flux tube 
refills with a cold plasma. The study was performed using simulations based on a small-scale PIC 
code. The results of this study are summarized in two companion papers which are submitted for 
publication in the Journal of Geophysical Research [Singh and Leung, 1994a, b]. In Paper 1 
[Singh and Leung, 1994a] we have studied the role of ion-beam driven instabilities in filling a 
magnetic flux tube. Large-scale models of plasmaspheric refilling have revealed that during the 
early stage of the refilling counterstreaming ion beams are a common feature. However, the 
instability of such ion beams and its effect on refilling remain unexplored. The difficulty with 
investigating the effect of ion-beam driven instability on refilling is that the instability and the 
associated processes are so small-scale that they cannot be resolved in large-scale models; 
typically the instabilities have scale lengths of a few tens of plasma Debye length, which is a few 
meters at the most, and the spatial resolution in large-scale models is at least several tens of 
kilometers. Correspondingly, the temporal scale of the instability is by several orders of 
magnitude smaller than the temporal resolution afforded by the models. In order to learn the 
basic effects of ion beam instabilities on refilling, we have performed numerical simulations of the 
refilling of an artificial magnetic flux tube. The shape and size of the tube are assumed so that the 
essential features of the refilling problem are kept in the simulation and at the same time the 
small-scale processes driven by the ion beams are sufficiently resolved. Two types of simulations 
have been performed, in one type we treat ion kinetically and electrons are assumed to obey the 
Boltzmann law. In the other type both electrons and ions are treated kinetically. A comparison 
between the results from such simulations reveal that in the latter type of simulations electron-ion 
(e-i) and ion-ion (i-i) instabilities occur and significantly modify the evolution of the plasma 
density distributions in the flux tube along with the total plasma content. When the electron 
dynamics is simplified by the assumption of the Boltzmann law, both the electron-ion and ion-ion 
instabilities are inhibited, and only in very late stage of the filling there is a weak scattering of ions 
due to an enhanced plasma fluctuation level. On the other hand, when electrons are treated 
kinetically, the e-i instability occurs at an early stage when ion beams are too fast to excite the i-i 
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instability. The former instability heats the electrons so that conditions for the latter instability are 
eventually met. The i-i instability and its non-linear evolution creates potential structures 
including several electrostatic shock pairs which significantly modify the filling process. The 
electrostatic potential structures are highly dynamic, and at times they appear as moving double 
layers greatly affecting the state of the plasma inside the central region of the flux tube. 

In Paper 2 [Singh and Leung, 1994b] we have studied the effects of equatorially trapped hot 
plasma on refilling. Equatorially trapped hot plasmas are a common feature of the outer 
plasmasphere, where flux tube refilling with cold ionospheric plasma occurs after magnetic 
storms. The role of the hot plasma consisting of hot anisotropic ions and isotropic warm 
electrons in the refilling process is examined by means of numerical simulations using a one- 
dimensional particle-in-cell code. Simulations are performed on the filling of an artificial flux tube 
having a minimum magnetic field at its center. We have performed two types of simulations; in 
one type, called here Run-A, we allowed cold plasmas to flow into a centrally trapped hot plasma 
consisting of warm isotropic electrons and hot anisotropic ions with perpendicular temperature 
T±>T\\, the parallel temperature. Run-A reveals a variety of plasma processes relevant to the 
plasmaspheric refilling affected by the presence of a hot plasma, including formation of 
propagating electrostatic shocks, intrinsically unstable plasma distribution functions produced by 
the mixing of hot and cold plasmas, weak downward electric fields supported by an extended 
potential distribution in the relatively late stage of the evolution of plasma in the flux tube, and an 
enhanced flux tube filling. In the other type of simulation first a cold plasma flow was allowed to 
set up in the flow, then a hot plasma consisting of the isotropic electrons and anisotropic ions 
(Tl > ?ii) was suddenly injected into the central region of the flux tube. In this case the main 
distinguishing feature was the formation of relatively stable shocks near the mirror points of the 
centrally trapped hot plasma. The shocks were found to be standing, unlike in the previous type 
of simulation. Since the standing shocks form near the effective mirror points of the centrally 
trapped hot ions, they are called mirror shocks to contrast them from the moving electrostatic 
shocks seen in Run-A. The stability of the standing shocks was found to increase with the 
decreasing temperature of the warm electrons injected with the hot plasma. Wherever possible, 
similarities between the results from the simulations and those from observational data are 
pointed. 

Tasks for the Grant Period Beginning August 1, 1994 : 

After the investigations based on small-scale simulations reported in Papers 1 and 2, we 
propose to include the processes seen in these simulations in large-scale hydrodynamic and 
semikinetic models for the plasmaspheric refilling. We will first attempt to include the effects of 
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electron-ion instability in terms of electron heating. For this purpose we must develop an 
algorithm for the electron heating rate as the ion beams evolve during the refilling. After this 
phase of the work, we will attempt to include the effects of ion-ion instability on the plasma. 
Since the ion-ion instability leads to the formation of vortices and electrostatic shocks, it is not 
clear at this time how to include its effect in the large-scale models. We propose to investigate 
this issue. 

One of the interesting results from the small-scale simulations is that the mixing of cold 
plasma with the equatorially trapped hot plasma produces low energy ring distributions in the 
perpendicular velocity. Such a distribution function is known to excite the lower hybrid and ion- 
Bemstein waves. We propose to investigate the stability of the distribution functions seen in the 
one-dimensional model by performing 2-dimensional PIC simulations, which will allow us to study 
transverse heating of the cold ions. Our main goal here will be to critically examine the conditions 
under which the ion heating produces the trapped ion population merged with a cold core as 
observed from DE-1 [Olsen et al, 1987], 
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PLASMA: ELECTROSTATIC SHOCK FORMATION , 

^3 A 38^7 

Nagendra Singh ' 
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Abstract . Effects of equatorially trapped hot plasma on 
the highly supersonic cold-plasma flow occurring during 
early stage plasmaspheric refilling are studied by means 
of numerical simulations. It is shown that the equatori- 
ally trapped hot ions set up a potential barrier for the 
cold ion beams and facilitate formation of electrostatic 
shocks by reflecting them from the equatorial region. Sim- 
ulations with and without the hot plasma show different 
flow properties; the formation of electrostatic shocks oc- 
cur only in the former case. The simulation with the hot 
plasma also reveals that the magnetic trapping in con- 
junction with the evolution of the electrostatic potential 
barrier produces ion velocity distribution functions con- 
sisting of a cold core and a hot ring in the perpendicular 
velocity. Such a distribution function provides a source 
of free energy for equatorial waves. The corresponding 
electron population is warm and field aligned. 

Introduction 

So far most theoretical studies on plasmaspheric refill- 
ing have been primarily concerned with the outflow of cold 
ionospheric plasma and its trapping in the flux tube. In 
such theoretical studies an important observational fact, 
which has been ignored, is that the flux tubes undergo- 
ing refilling contain a hot plasma population trapped in 
the equatorial region. Such plasmas originate from the 
ring current or the plasma sheet and are characterized by 
TK > T$, where Tj[ and I’jf are the perpendicular and 
parallel temperatures of the hot ions, respectively. Spe- 
cific observations of such hot ions having relatively large 
pitch angle anisotropies (A t = T*^/T^ > 2) come from 
G EOS-1 and -2, which observed the hot ions in the en- 
ergy range > 10 keV in the noon sector. Such hot ions 
are known to excite electromagnetic ion cyclotron waves [ 
Roux et al., 1982]. 

Another set of observations, where the role of hot 
plasma has been invoked, deals with thermal ions trans- 
versely heated to energies up to a few hundred eV and 
trapped in the equatorial region [Olsen et al, 1987]. Here 
again hot ions are suggested to be the source of free en- 
ergy for exciting the broadband waves observed from DE- 
1. Theories suggest that the broadband waves are driven 
by a combined effect of temperature anisotropy and ring 
type of distributions of energetic protons [Perraut et al, 
1982 ]. The effect of heating of the thermal ions on re- 
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filling has been studied [ Singh and Chan, 1992], but the 
direct effect of the hot ions on the refilling has not been 
studied so far. 

The purpose of this letter is to show that the hot 
anisotropic ions, which are commonly present in the equa- 
torial region of flux tubes undergoing refilling and drive 
the equatorial processes discussed above [ Olsen 
et al, 1987; Roux et al, 1982], facilitate the process of 
electrostatic shock formation since they provide an effec- 
tive potential barrier for the upfiowing ion beams of cold 
ionospheric plasma. The shock formation occurs even for 
highly supersonic ion beams which are expected to occur 
during early stage refilling [Banks et al, 1971; Singh et al, 
1986; Rasmussen and Schunk, 1988; Singh , 1990; Wilson 
et al, 1992]. If the hot plasma consists of isotropic elec- 
trons and anisotropic biMaxwellian ions with anisotropy 
Ai = T*/T$ > 1, the potential difference between the 
equator, where the minimum magnetic field strength is 
B m , and an ionospheric point where the magnetic field 
is B , is given by $ = kT^ / e[ 1 -f T t ^/T e ]~ l ln(T), where 
T = Ai(l - B m /B) + B m /B [ Whipple , 1977]. 

Since B m /B << 1, the equatorial potential with respect 
to the ionosphere for T e << is $ ~ (kT e /e) ln( A t ). For 
A, = 2 and kT e /e = 10 V, the typical values for these pa- 
rameters used in wave analysis [Roux et al, 1982], $ ^ 
7 V. However, this potential difference is true when no 
ionospheric cold plasma is present; in the presence of a 
cold plasma, it is not certain how large the potential dif- 
ference is and how it is distributed along the field line. 
In the following discussion, we present results from nu- 
merical simulations elucidating the interactions between 
cold and hot plasmas when the former plasma flows into 
the latter one trapped in a magnetic mirror. The simula- 
tion presented here is small scale, in contrast to the large- 
scale problem in space. Therefore, the results presented 
here only serve the purpose of elucidating the processes 
involved in hot-cold plasma interactions. 

Numerical Model 

We perform a one-dimensional particle-in-cell simula- 
tion of plasma flow along an artificial flux tube (Figure 
la). The magnetic field B(X) = B 0 (l - a exp[— (A — 
d/2) 7 /a 2 )) where B 0 is a constant field outside 
the minimum-field region, d is the size of the simulation 
system and the choice of a and a determines the desired 
field distribution. The hot plasma is created by injecting a 
large number of electrons and ions in the minimum-B field 
region. Such plasma particles are chosen from Maxwellian 
or bi-Maxweilian velocity distributions depending upon 
the requirement of a simulation run. The cold plasma 
flows into the flux tube from the two plasma reservoirs at 
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the ends of the simulation system at A" = 0 and X = d 
(Figure la). The simulation technique is described in 
Singh and Chan [1992]. We solve for the motions 
of charged particles in self-consistent electric fields, de- 
termined by solving Poisson’s equation. As the particles 
move in the flux tube, their magnetic moments are as- 
sumed to be conserved. 

In the cold-plasma reservoirs electron and ion temper- 
atures are TV The reservoirs supply a continuous flux of 
charged particles into the flux tube through the process 
of plasma expansion. In the simulations reported here, 
we have used m.,- / T 7 i e = 400, which adequately separates 
the electron and ion time scales and, at the same time, 
allows computationally feasible runs. The hot plasma is 
injected into the central region of the simulation system 
at time t=0 (Figure la). The properties of the hot plasma 
are described later. Numerical parameters of the simula- 
tions are as follows: system size d = 5 xlO 3 B field 
parameters a = 0.9, cr — 750 A^, cell size Ax — 20 Aj, and 
time step At = O.lu;^, where A<* and are the Debye 
length and electron-plasma frequency in the cold plasma 
reservoirs. In the following discussion we have used nor- 
malized quantities defined as follows: time f = tu)^, dis- 
tance X = X /\ <*, velocity V = VjV te and electric poten- 
tial $ = e$/kT 0 , where V it — (kT Q jm e ) 1/2 and Wpi 0 = 
(m e /mj 

Numerical Results 

First we present results from a simulation in which 
“equatorial” hot plasma is not included. This simula- 
tion serves as a reference against which the hot plasma 
effects on the cold plasma flow are compared. Figures lb 
to If show the temporal evolution of the flow of cold ions 
into the flux tube from the reservoirs shown in Figure la. 
These panels show the phase-space plots in the X - Vjj 
plane. Each dot in the figure represents an ion, giving 
its parallel velocity and position. The simulation begins 
at t = 0, when plasmas from the reservoirs begin to flow 



x(1000A*) 

Fig. 1. (a) Artificial flux tube in which plasma expands 
from the cold plasma reservoirs at X=0 and X d. The 
hot plasma is injected in the central region, (b) to (f) 
ion phase-space plots in X - Vy plane at selected times as 
shown. The plots show evolution of the cold plasma flow 
in absence of a hot plasma. 


in the flux tube. The corresponding flow of electrons is 
not shown here. Figure lb shows expanding ion beams 
into the simulation region at f = 200. The plot for t = 
1,000 (Figure lc) shows that the ion beams have expanded 
from one end to the other end of the simulation region, 
setting up counterstreaming. Equatorial shock formation 
has not occurred. At this time, ion beams are too fast 
to couple through the ion-ion instability; the relative flow 
velocity between the beams is V re i = 0.4 V* — 8c # , where 
c, is the local ion-acoustic speed at the equator, X = d/2, 
and is about .051*. However, at later times the beams 
become progressively slower as the plasma accumulates in 
the flux tube, and when they become sufficiently slow in 
the regions near the ends of the flux tube they excite ion- 
ion instability as seen for l > 1,400 (Figures Id to If); 
the instability creates vortex-like structures in the phase 
space. 

The instability is seen to thermalize the ion beams (Fig- 
ure le). However, in the middle of the simulation region, 
where the magnetic field is minimum, counterstreaming 
of ion beams is seen to persist. 

Simulation with a hot plasma in the region of mini- 
mum magnetic field reveals a quite different behavior of 
the flow. We performed several simulations by varying the 
temperatures of the hot plasma but kept ion anisotropy 
Ai >1 and electron anisotropy A t — 1 and T t << T^. 
Such properties of hot plasma have been measured [Roux 
ct al. t 1982]. Relative densities of the cold (nc) and hot 
(n/f) ions in their respective source regions were also var- 
ied. The results from such a parametric study will be re- 
ported elsewhere. We present results here for — Tf^/2 
= 900To, T t = lOTo, and n c — 100 

The hot plasma is injected in the central region 2300 < 
x/\ d < 2700 at t = 0, when cold plasma flows begin from 
the reservoirs. Figure 2 shows the evolution of the cold 
plasma flow along with the potential and density profiles 
in the flux tube. The panels a to d show the ion phase 
space plot in X — Vjj plane. The phase-space plot of ions in 
X — Vx plane is shown in panels e to h. The evolution of 
the potential profiles in the flux tube is shown in panels i to 
1, and the corresponding plasma density is shown in panels 
m to p. The evolution of the electron phase- space (X - Vj|) 
is shown in panels q to t. As shown on the top of the figure, 
the columns from left to right show the properties of the 
flow at i - 100, 500, 700, and 1000, respectively. At an 
early time (t = 100), the flow can be characterized by the 
following features: (1) expanding cold ion beams (panel 
a), (2) equatorially trapped ions originating from the hot 
plasma (panel e), (3) a potential maximum at the mid- 
point of the simulation region (panel i), and (4) a density 
maximum coinciding with the potential maximum (panel 
m). The density maximum at this time entirely represents 
the trapped hot ions, and there is no contribution from the 
cold plasma flow. 

At t — 500, the cold plasma flow comes into contact 
with the hot plasma; the cold ion beams have somewhat 
penetrated into the region of the hot plasma (panels b 
and f). Already the cold beams show the sign of retar- 
dation by the potential barrier. The slowing down of the 
ion beams enhances the potential barrier further; panel j 
shows that at t - 500 the equatorial potential maximum 
has grown to 15fcTo/e, in contrast to 3fcT 0 /e at t = 100. 
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Fig. 2. Evolution of the cold plasma flow in the presence of equatorially trapped hot ions. Each column shows the 
state of the plasma at times shown at the top. (a) to (d) X - Vjj phase-space plots for ions, (e) to (h) X - Vi phase 
space plots for ion, (i) to (l) Potential distributions, (m) to (p) Ion density distributions, and (q) to (t) Electron 
phase-space plots in X — V[j plane. 


Note the change in the vertical scale between panels i and 
j. The filling of the flux-tube with the cold plasma changes 
the density profile by creating a density minimum at the 
equator (panel n). 

The plots at t = 700 show that the contact between 
the hot and cold plasma has evolved into a pair of elec- 
trostatic shocks, one on each side of the hot ions trapped 
in the minimum-B field region. The shocks are indicated 
by arrows where flow velocities of the cold ion beams sud- 
denly decrease (panel c). Near the shocks, the potential 
profile shows sharp jumps (panel k) and the density profile 
shows spikes (panel o). 

The shocks propagate away from the “equatorially” 
trapped hot plasma towards the end of the flux tube, 
thermalizing the cold plasma just behind the shocks. The 
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Fig. 3. Ion velocity distribution functions in the equato- 
rial region (a) Parallel velocity distribution, (b) Perpen- 
dicular velocity distribution. 




propagation of shocks can be seen by comparing the pan- 
els in columns for i = 700 and 1000. The shocks propa- 
gate with a speed V,h 2 ^ 0.09 2 ; C t4>> where C to is the 

ion-acoustic speed in the cold plasma reservoirs. As the 
shocks propagate away from the hot plasma, the thermal- 
ized cold plasma behind the shocks punches through the 
region of hot plasma where the potential barrier existed 
earlier. This penetration of the cold plasma is clearly 
seen from the X - Vj_ plots shown in Figure 2e to 2h. 
As the cold ions penetrate into the minimum-B field re- 
gion, where hot ions are trapped, they cool adiabaticlly 
as clearly seen from these plots. The penetration of cold 
plasma is accompanied with a reduction in the potential 
barrier, caused by the cold electrons being accelerated into 
the high potential region (panels r and s); these electrons 
tend to neutralize the effect of the hot ions stably trapped 
in the minimum-B field region. 

A noteworthy feature of the plasma after the spatial 
mixing of hot and cold plasmas is the nature of the re- 
sulting ion distribution functions in the equatorial region, 
as shown in Figures 3a and 3b. The parallel velocity dis- 
tribution function (Figure 3a) shows a core of cold ions 
superimposed on a warm ion population. The perpendic- 
ular velocity distribution function (Figure 3b) also shows 
the cold ions, but the hot ions appear as a beam. Since 
the ions in the beam are nearly uniformly distributed in 
their phase, the beam is actually a ring in the perpendic- 
ular velocity space. The ring distribution function is the 
result of the trapping of the hot ions in the magnetic mir- 
ror in combination with the evolution of the electrostatic 
potential distribution. Some of the ions from the hot pop- 
ulation having relatively small perpendicular velocities are 
lost from the mirror due to the parallel electric fields. In 
the absence of the electrostatic potential, one expects a 
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loss-cone distribution. The loss of ions with relatively low 
VI, yields the ring distribution. The cold core originates 
from the cold ions penetrating in the minimum- B field re- 
gion. The importance of these findings lies in the fact that 
ring types of distribution functions along with a core of 
cold ions can be a source of free energy for exciting equa- 
torial waves [Lee and Birdsall, 1979], which can heat the 
latter ions. We point out that the low energy end of the 
ring distribution found here occurs at an energy of a few 
tens of eV, unlike the relatively energetic ring distribution 
with energy > 5 keV observed from GEOS [Perraut et ai, 
1982]. 

Conclusion and Discussion 

The main conclusions of the study in this paper are the 
following: (1) The presence of the hot anisotropic ions 
trapped in the equatorial region can drastically affect the 
flow of the cold plasma. (2) The potential barrier associ- 
ated with such a hot ion population facilitates formation 
of electrostatic shocks even when the cold ion beams com- 
ing from the ionosphere are highly supersonic. (3) The 
eventual mixing of the hot and cold plasmas produces a 
potentially unstable velocity distribution function for the 
ions. They develop a ring distribution in the perpendic- 
ular velocity. The ring distribution with df(V±)/dV± be- 
gins at a relatively low energy corresponding to the energy 
of an ion falling down the potential barrier of the order of 
10 V. The ring distribution discussed here is different from 
the observed high energy rings (> 5 keV) produced by hot 
plasma injection and its subsequent convection [ Perraut et 
ai, 1982]. 

Recent satellite observations have shown that equato- 
rially trapped warm ions are generated by equatorial per- 
pendicular heating of the thermal ions; the heating is 
caused by waves generated by an energetic hot ion popula- 
tion, having the temperature anisotropy and/or 

ring type of velocity distribution [ Olsen et al ., 1987; Per- 
raut et ai, 1982 ]. Olsen et ai [1987] also report the 
reflection of cold ion beams. However, it is believed that 
the reflection is caused by the potential barrier set up 
by the thermal ions which have undergone the transverse 
heating. From our simulations, it appears that the hot 
ion population, which drives the waves by virtue of its 
anisotropy or ring, may reflect the fast ion beams at an 
earlier stage before the waves grow and consequent ion 
heating produces a warm anisotropic ion population. 

Futhermore, before thermal ions are heated they must 
enter the hot plasma region. This entry is retarded by the 
potential barrier set up by the hot ions. Even if the hot 
plasma injection occurs at a stage when field-aligned flows 
have set up at the equator, the potential barrier is likely 
to expel the cold plasma from the equatorial region. The 
presence of thermal ions is essential in the wave generation 
processes. Therefore, it is suggested that the reflection of 
ion beams observed from DE-1 may not be necessarily 
caused by the equatorial perpendicular heating of ther- 
mal ions; it appears that the reflection of ion beams and 
heating of thermal ions are all driven by the hot ions, but 


heating should phenomenologically follow the reflection of 
upflowing ion beams. 
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Abstract. Hydrodynamic and semikinetic treatments of plasma flow along closed geomagnetic 
field lines are compared. The hydrodynamic treatment is based on a simplified 16-moment set 
of transport equations as the equations for the heat flows are not solved; the hear flows are 
treated heuristically. The semikinetic treatment is based on a particle code. The comparison 
deals with the distributions of the plasma density, flow velocity, and parallel and perpendicular 
temperatures as obtained from the two treatments during the various stages of the flow subject 
to certain assumed boundary conditions. In the kinetic treatment, the appropriate boundary 
condition is the prescription of the velocity distribution functions for the particles entering the 
flux tubes at the ionospheric boundaries; those particles leaving the system are determined by 
the processes occurring in the flux tube. The prescribed distributions are half-Maxweiiian with 
temperature T 0 and density n Q . In the hydrodynamic model, the prescribed boundary 
conditions are placed on density (« 0 ), flow velocity (v o ), and temperature (T 0 ). We found that 
results from the hydrodynamic treatment critically depend on v o ; for early stages of the flow 
this treatment yields results in good agreement with those from the kinetic treatment, when 
v o » (kT 0 / 2 xm) U2 , which is the average velocity of particles moving in a given direction for a 
Maxwellian distribution. During this early stage, the flows developing from the conjugate 
ionospheres show some distinct transitions. For the first hour or so, the flows are highly 
supersonic and penetrate deep into the opposite hemispheres, and both hydrodynamics and 
kinetic treatments yield almost similar features. It is found that during this period heat flow 
effects are negligibly small. When a flow penetrates deep into the opposite hemisphere, the 
kinetic treatment predicts reflection and setting up of counterstreaming. In contrast, the 
hydrodynamic treatment yields a shock in the flow. The reasons for this difference in the two 
treatments is discussed, showing that in view of the relatively warm ions, the coupling of ion 
beams and the consequent shock formation in the off-equatorial region are not likely due to the 
enhancements in the beam temperatures. The counterstreaming in the kinetic treatment and 
the shock in the hydrodynamic treatment first advance upward to the equator and then 
downward to the ionospheric boundary from where the flow originated. The transit time for 
this advancement is found to be about 1 hour for the perspective models. After 2 hours or so, 
both models predict that the flows from the ionospheric boundaries are generally subsonic with 
respect to the local ion-sound speed. At late stages of the flow, when a substantial fraction of 
ions entering the flux tube begin to return back in the kinetic treatment, the hydrodynamic 
treatment with the boundary condition v o = ( kT Q I2xm) ui yields an overrefilling, and the 
choice of v o becomes uncertain. 


1. Introduction 

In connection with the problem of plasmaspheric refilling, in 

recent years several models for plasma flow along closed 

magnetic field liens have been developed [Khazanov et a/., 

1984; Singh et al f 1986; Singh et al, 1988; Rasmussen and 
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Schunk , 1988; Singh et aL, 1991; Gutter and Gombosi , 1990; 
Wikow ef aL 1992]. These models differ in complexity in 
tpr m * of describing the plasma and in including the ionosphere 
as a source of plasma for the refilling. In terms of describing the 
plasma »nd the most contrasting feature of the existing models 
with the hydrodynamic and kinetic treatments for the 
flows, based on plasma fluid equations and a particle- in-cell 
code, respectively. For the purpose of including the ionosphere 
as a source of plasma for the refilling, in most studies the 
topside ionosphere is r ep laced by a set of boundary conditions 
on the plasma flow, except by Guiter and Gombosi [1990], who 
have included the generation and loss of plasma through 
chemical reactions in a hydrodynamic model. In this paper, we 
are mainly concerned with the treatments of the plasma with 
simple sets of boundary conditions on the plasma flow, we 
compare the properties of the plasma along closed field lines as 
given by hydrodynamic [Singh et a/, 1992] and semikinetic 
[Wilson etal , 1992] treat m e nt s. 

The s uccess of a hydrodynamic treatment depends on the 
problem being solved and on the ingenuity of the user in 
choosing the hierarchy of moment equations on which fluid 
equations are based. In recent years, researchers in space 
physics have used fluid descriptions based on 1 3-moment 
[Schunk, 1977; Mitchell and Palmadesso , 1983], and 16- 
moment [Barakat and Schunk ; 1982; Ganguii and Palmadesso , 
1987; Gombosi and Rasmussen, 1991; Korosmerzey et a/., 1992, 
1993] set of transport equations. In developing the moment 
equations, the ingenuity Lies in a series expansion of the plasma 
distribution f unc tion using a biMaxweliian distribution function 
as a base. Therefore hydrodynamic treatment based on moment 
equations are good as long as the distribution function is close to 
a biMaxweliian. When the distribution function severely departs 
from a biMaxweliian and involves multistreaming of plasma 
particles, the moment equations are seriously handicapped, 
despite the sophistication of the higher-order moment equations 
used. 

As mentioned above, recent models for plasmaspheric 
re filling are based on both a kinetic treatment using particle-in- 
cell (PIC) code and a hydrodynamic treatment with varying 
degrees of sophistication in choosing the hierarchy of the 
moment equations. In some early models, only continuity and 
momentum equations were solved for the ions, and the electrons 
were flggnrrayl to remain isothermal [Singh et aL, 1986; 
Rasmussen and Schunk , 1988; Singh, 1988]. Studies that 
included temperature equations assumed that either the heat 
flow is given by the colhrion-ckxninated thermal conductivity 
[Khazanav et aL, 1984; Guiter and Gombosi, 1990] or ignored 
the heat flow completely [Singh, 1992; Singh and Chan , 1992]. 
Neither of these assumptions correctly describe the heat 
transport in the refilling problem and Horwitz , 1992]. In 
the coilisionless limit of plasma flow during refilli ng , the usual 
description of heat flow in terms of Spitzer thermal conductivity 
breaks down, and such a treatment overestimates the heat flow. 
When the heat flow becomes large, the validity of such 
equations ceaw«- and numerical instabilities result in 
computational work. Since , a priori it is not known when a 
large heat flow develops in a model treating plasma flow in a 
flux tube extending to altitudes of several Earth radii, ad hoc 
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damping m echani sms are included to damp the heat flow 
whether it is physically warranted or not [Palmadesso et ai, 
1988; also C. E. Rasmussen, private communications, 1993]. In 
modeling of plasma flow at relatively low altitudes, where 
coilisional effects are important, the above problem with the 
heart flow does not appear to be of a major concern [Korosmezey 
etal, 1992]. 

Despite the above difficulties with the hydrodynamic 
treatment, it has been used for practical reasons because it 
provides simplicity and considerable economy in computational 
work, and depending on the plasma conditions it can work 
successfully. Therefore it is advisable to keep in mind the 
assumptions made in using this treatment and, if possible, it is 
even better to check the validity of this treatment by comparing 
its prediction a gainst that from a kinetic treatment Such a 
comparison may reveal when and how a fluid model succeeds. 

The purpose of this paper is to cany out a comparison 
between models of the plasma flows along closed field lines 
based on kinetic and hydrodynamic treatments. The former 
treatment uses a PIC code for ions [Wilson et a/., 1992], The 
latter one uses transport equations for the flow of mass, 
momentum, and parallel and perpendicular temperatures of ions 
[Singh, 1992], but the heat flow is treated heurisdcally { Metzier 
et aL , 1979]. In both the treatments, electrons are assumed to 
obey the Boltzmann law. In the present paper the ionospheric 
outflows is included by imposing a set of boundary conditions on 
the flow of ions at an altitude of 2000 km. The choice of this 
altitude is primarily due to the existing models [Wilson et al. 9 
1992; Singh, 1992] in which ionospheric loss and generation 
processes for the plasma are not yet included. 

The closed field lines provide the possibility of a variety of 
flow conditions ranging from highly supersonic to subsonic 
flows as an empty flux tube refills. Furthermore, the flows along 
closed fields lines develop counterstreaming due to 
interhemispheric plasma flows. Since hydrodynamic treatments 
are most suspect under counterstreaming situations [Manheimer 
et a/., 1976], the comparison earned out here provides a useful 
guide for assessing the validity and usefulness of a 
hydrodynamic treatment 

We have found that for the conditions of highly supersonic 
flows, the two-stream hydrodynamic treatment yields flow 
properties in good agreement with that from the senukinetic 
treatment Demars and Schunk [1991] reported a similar 
agreement based on 16-moment set of equations including heat 
flows. We find that the bulk parameters such as the density, 
flow velocity, and temperatures are in good agreement even for 
the simpler hydrodynamic model when heat flow is included 
heurisdcally; the reason being simply that when the flow is 
highly supersonic, the dominant transport of heat is through the 
bulk flow velocity and the transport due to the thermal effects is 
negligibly small. 

When reflection of flows causes counterstreaming, the 
hydrodynamic treatment gives rise to shock formation, which is 
not seen from the kinetic treatment An examination of the 
shock formation through ion-ion instability shows that off- 
equatorial shock formation is inhibited by an unfavorable 
temperature conditions on electron to ion temperature ratio. The 
issue of the equatorial shock formation remains unsettled here 
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due to the coa rseness of the spatial and temporal resolutions in 
both the kinetic and hydrodynamic models. When the 
counterstreaming flows become subsonic, the hydrodynamic and 
semikinetic treatment again produce flow properties in 
reasonable agreement 

In order to study plasma flow in space, the plasma treatment 
must be supplemented by a set of boundary conditions on the 
flow equation. For the closed field lines, the boundary 
conditions are determined by the topside ionosphere. The 
boundary conditions along with the demand for plasma at high 
altitudes produce the flow. The ionospheric boundary conditions 
involves generation and loss of ionospheric plasma particle 
species. Since here our primary goal is in identifying the kinetic 
and fluidlike behaviors of plasma flow and not the supply of 
plasma from the ionosphere and the refilling rate, we have 
simulated the outflow of ionospheric plasma by imposing a set 
of boundary conditions at an altitude of 2000 km in both the 
hemispheres. In the semikinetic model, the imposed boundary 
condition is on the velocity distribution function of the ions 
entering the flux tube. It is assumed to be half-Maxweilian. The 
returning particles are self-consistently determined. In the 
hydrodynamic model the boundary conditions are the moments 
of such a distribution. Since the kinetic effect dealing with the 
returning particles are lost in the hydrodynamic model, the 
hydrodynamic model does not agree with the kinetic model 
when returnin g ion flux becomes sufficiently large. Can this 
disagreement be resolved by a more sophisticated treatment of 
the heat flow by using a complete set of 16-moment equations? 
In order to answer this question, further comparative studies are 
suggested. 

The rest of the paper is organized as follows. The theoretical 
models are described in section 2. The comparison between the 
results fr om the two models is carried out in section 3. The m a in 
conclusions of the paper and their discussion are given in 
section 4. 


2. Theoretical Models 

The semikinetic model, which is based on a particle-in-cell 
code, has been previously described for both open [Wilson et a/., 
1990] and closed [Wilson et ai , 1992] flux tubes. Coulomb 
collisions are included in the model, the collisions are 
implemented by pairing simulation ions according to an 
algorithm which conserves energy and momentum [Takizuka 
and Abe , 1977]. The algorithm yields good approximation for 
the collisions when the coilisional relaxation tune is shorter t han 
the time step in advancing the ion motion. In the hydrodynamic 
model, we solve the plasma transport equations based on 16- 
moment approximation [e.g., Ganguli and Palmadesso , 1987, 
and Barakat and Schunk, 1 982] : 


dn d 1 dA 

— +— (nV) = -nV— — 
dt ds A ds 


( 1 ) 


3P d e 3T\ { 1 dn 

dt ds m ds n ds 


- 4 - 



1 dA & 

-X»(r) ~(k/m)(T t -T 1 )——+ [—], 
A ds at 


( 2 ) 


at 


^ [VT ' ] = ~ T * 


£1 

ds 


i i a 

n A as 


2 1 0A 

{qnA)+—q L — — + 
n A as 



£L 

at 


a av 

+— ^rj-r A — -r x r- 


1 

7 ^ 


i 

fl 


1 A4 
/I <2* 


(3) 


i i a 

n A 0s 


<*i ^)+ 



(4) 


where t is time, r is the geocentric distance along the flux 
tube, j is the distance along the tube from the northern 
ionospheric boundary at X = +A o (see Figure 1); /r, K, r H , and 
r x are the number density, flow velocity, and parallel and 
perpendicular temperatures of ions in the plasma flow, 
respectively ; q t and q L are the heat fluxes along the magnetic 
field line associated with T n and T lt respectively; £ is the 
parallel electric field; g H is the component of the gravitational 
force parallel to the magnetic field, and m and e are the ion 
mass and charge, respectively. The collision terms denoted by 
[]* are calculated using Burger's Formulae [Burger, 1969], 
which are modified to include flow velocity corrections [ Mitchell 
and Palmadesso , 1 983; Ganguli and Palmadesso , 1 987] and the 
correction for temperature anisotropy [Ichimam et al. , 1973; 
Singh, 1991]. 

We do not solve the heat flow equations, which have proven 
to be quite troublesome to solve numerically [Palmadesso et al 
1 988]. The difficulty arises for a relatively large heat flow, for 
which the moment equations themselves become invalid. Since 
it is unpredictable in a model when the heat flow may be large, 
ad hoc procedures are employed to attenuate the heat flow for 
the numerical stability of the models. This has been found to be 
true irrespective of the numerical techniques used for solving the 
equations [Palmadesso et al 1988; Korosmezey et al. , 1992, 
1993; also C. E. Rasmussen , private communication, 1993]. 

In this paper, instead we have included the effects of heat 
flow heunsticaily by closely following the treatments in solar 
wind studies [e g., Metzler et al. , 1979]. In a collisionless 
plasma the usual picture of heat flow, given by q m - - K a VT a 
with K a as the thermal conductivity, may not be valid because 
L, the mean free path, is » L? = (T~ l cT f ds)' 1 , the scale 
length in the temperature variation. In such a collisionless 
situation, the heat flux can be calculated on physical ground as 
follows. The unidirectional heat fluxes and across a 
surface in a plasma described by a biMaxweilian distribution 
function with parallel and perpendicular temperatures T n and 
T l , respectively, say along the magnetic field vector, are given 
by 


<?« =nkT^kTil2xmY n (5) 

ql =nkT l (kT t /2xmy n (6) 
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In a unif orm plasma for which the distribution function is 
independent of the parallel coordinate, the heat flow at any point 
is zero because he at flux in a given direction is c ance lled by the 
heat flux in the opposite direction. In the presence of a spatial 
inhomogeneity, the cancellation is not complete, and the heat 
fluxes q n and q x appearing in equations (3) and (4) an be 
heuristically written as {Metzler ex aL, 1 979] 

q a =stinkT a V, t (7) 

where the subscript a stands for 1 or II , s= -1 if 
dT a !da> 0, and s= 1 if JT a /ds<0 . Thus in the heat flow 
model adopted here, only the sign of the heat flux depends on 
the temperature gradient and not its magnitude. The factor 77 
determines the reduction in the heat flow below the 
unidirectional fluxes in (5) and (6). Later on in this paper we 
show that t) in the range say 0.1 - 0.3 yields results in a 
reasonable agreement with the kinetic model, in which heat 
fluxes appear self-consistently. A similar model for heat flow 
was used by Singh [1992] for plasma flow along open field lines, 
hi both the hydrodynamic and kinetic models adopted here, 
electric field E is calculated by assuming that the electrons obey 
the Boltzmann law and the condition of quasi-neutrality 
prevails. 

The plasma flow along a closed Geld line is studied by 
solving and initial -boundary value problem. In the hydrodynamic 
model, the plasma flows originating from the conjugate 
ionospheres are treated as separate fluids; this t rea t m e nt is 
termed as a two-stream model [Smg/i, 1988; Rasmussen and 
Schunk ; 1988; Singh , 1990]. In both the models, it is assumed 
that at the initial time (f = 0) the flux tube is highly depleted. 
The depletion is given by n { =« c (sinA/ sinA 0 ) 6 , with the 
minimum density limited to 10^n o , where n 0 is the density at 
the ionospheric base X = ±X Q (Figure 1). Initial flow velocity 
V(X t f = 0) = 0 and temperature T L (Xj = 0) = T^Aj = 0) = 
T 0 = 0.3 eV . In the hydrodynamic model, the boundary condition 
for the fluids originating from the northern hemisphere are 
n m {k-JL 0 . t) = n ot V n (A = A gf t) = V ot W o ,t) = T,(A o J)=T 0t 
at the boundary X = -X 0 floating boundary conditions are 
applied. A set of similar boundary conditions is used for the 
fluid originating fro m the southern hemisphere, but with the 
roles of X » ±A 0 interchanged. In the kinetic model the 
boundary conditions on ion distribution function /(A, V) are 
that f(A = A 0 ,VZ 0) and f(X = -A oy V ZQ) are half- 
Maxwellians, with a temperature T 0 . These boundary conditions 
prescribe only the ions entering the flux tubes. The ions leaving 
the flux tubes are determined by the processes occurring inside 
it A half-Maxwell ian, and not a displaced Maxwellian, is 
chosen because of the following reasons. There is no clear 
observational evidence of supersonic flows along closed field 
lines at an altitude of 2000 km. Furthermore, our calculations 
show that in about 2 hours the flow in the flux tube becomes 
subsonic nearly all along the flux tube; only for an initial stage 
of about 2 hours, supersonic flows are seen. In view of such 
uncertainties on the flow velocity at 2000 km, a half-Maxwell ian 
serves the purpose of the comparative study. 
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3* Numerical results 


Wc compare here the properties of the flow in a flux tube 
with L = 4 as seen from the semikinetic and hydrodynamic 
models. Since the boundary value of the flow velocity (V o ) and 
the heat flow factor 7 in the hydrodynamic model are free 
parameters, comparison is performed by varying them over 
physically reasonable ranges. The comparison also deals with 
the accumulation of plasma in the flux tube and the equatorial 
plasma density. 

3.1 Initial Supersonic Flow 

We recall that the hydrodynamic model is based on two- 
stream flow in which flows originating form the two 
hemispheres are treated as separate fluids, and the temporal 
evolution of the two streams is separately studied. Likewise, 
even in the semikinetic model, the separate identity of the ions 
originating from the two hemispheres is maintained. Figure 2 
shows the evolution of the flow originating from the northern 
hemisphere as seen from the semikinetic model. This figure 
gives the phase-space density plots in X - V n plane, where X is 
the geomagnetic latitude and is the flow velocity along the 
magnetic field line. The positive and negative values of X 
correspond to the northern and southern hemispheres, 
respectively. The darkest region in the grey-scale plots 
represents the highest density as indicated by the scale on the 
right-hand side. At time t =0.003 hour, the plasma in the tube is 
essentially the initial plasma with a density profile given by 
n = nJsinA/ sinAJ 6 . At later times this plasma expands into 
the flux tube and it is seen to cross the equator at / =025 hour. 
Along with the expansion, new plasma enters the flux tube at 
the boundary X = X g . It is seen that by the time t = 0.75 hour, 
the flow has penetrated all the way to the opposite boundary at 
X = -X 0 . It is found that the plasma reaching this boundary is 
not totally lost, but it is partially reflected back, setting up a 
counterstreaming flow as seen from the plots for t ^ 1 hour. The 
reflected flow is seen to reach the boundary at X = X Q by the 
time t = 2 hours. The plasma flow originating from the southern 
hemisphere shows a similar behavior as shown in Figure 2, with 
the role of boundaries at X = ±X 0 interchanged. It is worth 
pointing out that after reflections, ions merge with the ion 
stream moving in the opposite direction, and they do not 
appear as a separate ion beam. 

The possible consequences of the counterstreaming flow will 
be discussed later on. We now compare the above features of the 
flow seen from the kinetic model with those seen from the 
hydrodynamic model. Figures 3a to 3d show the comparison for 
/ =0.5 hour, these figures show the distributions of (a) density, 
(b) flow velocity, (c) parallel temperature, and (d) perpendicular 
temperature. In each panel the curve from the kinetic model is 
labeled, and the curves from the hydrodynamic model for three 
values of the heat flow reduction factor 7 = 0, 7 = 0.05, and 
7 = 0.3 are indicated by the legend. It is seen from these figures 
that for most of the flux tube all four curves are quite close 
together, irrespective of the heat flow factor 7 . However, near 
the opposite boundary X = -A 0 , the curves from the 
hydrodynamic model tend to diverge from the kinetic model, 
depending on the value of 7 . This difference between the two 
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models is attributable to the feet when the flow begins to slow 
down due to a relative increase in the plasma density, the 
hydrodynamic model predicts an increase in the parallel 
temperature as clearly is seen for A < —36° in Figure 3c. The 
increase in T n enhances the pressure and further slows the flow 
and enhances the density. In the semikinetic model, the 
te mpe r a ture enhancement does not occur. Instead, a 
counterstreaming develops. This contrast between the two 
models becomes much clearer at later times, for example, at 
t - 1 hour for which the comparison is shown in Figures 4a to 
4d. 

Figure 2 shows that at / = 1 hour the reflected ions set up 
counterstreaming throughout the southern hemi sph ere (A < 0 ). 
Since the hydrodynamic model cannot handle the 
counterstreaming, the reflection process crea t es a shock, which 
is clearly ye n in the density, velocity, and temperature plots in 
Figures 4a to 4c, respectively, across the shock indicated by the 
arrows, density suddenly increases, flow velocity decreases and 
the parallel temperature also increases. We note that the 
hydrody namic curves for different values of tj begin to show 
some difference among themselves, with the curve for tj = 0.3 
being closest to that from the kinetic model. It is worth pointing 
out that the flow velocity in the kinetic model is the average 
over the counterstreaming ions. The average velocity is lower 
than that fro m the hydrodynamic model over the region of the 
counterstreaming, but where the counterstre aming has not yet 
occurred (A > 30*) the flow velocities from the two models are 
generally in good agreement 

The shock first propagates upward to the equator and then 
downward «nd reaches the ionospheric boundary at A = A 0 at 
t s 2 hours. The propagation of the shock in the density profile 
is shown by the arrows in Figure 5. The transit time of about 2 
hours for the shock is in agreement with the development of the 
rnimtpr gtr n*aming starting in the southern hemisphere and 
spreading to the northern ionospheric boundary by t = 2 hours, 
(see Figure 2). We find that the heat flow plays only a minor 
role in the motion of the shock; the shock speed is slightly 
pnhftnrgd with increased heat flow r \ ; for r] = 0.3 and 0.05 the 
shock is already absorbed near the boundary A = A 0 , while for 
r\ = 0, the shock can be still seen in the flux tube at t = 2 hours. 

After the shock reaches the boundary at A = A 0 , the flow in 
the flux tube becomes generally subsonic with respect to the 
ionacoustic speed, which is about 10 km / s -1 with electron and 
ion temperatures T 0 = 0,3eV at A = ±A 0 . We will discuss the 
subsonic stage after we examine the reason why a shock did not 
form during the early stage (~1 hour) of the counterstreazning 
(Figure 2) in the semikmetic model. 

3.1 Electrostatic Shock 

We have just seen that a shock automatically forms in the 
hydrodynamic model as soon as the flow begins to reflect near 
the opposite boundary. On the other hand, the kinetic model 
does not show the shock formation. Instead, a counterstreaming 
flow develops. We examine this issue in terms of the conditions 
for shock formation and ion velocity distribution function. 

According to the original suggestion of Banks et ai , [1971], a 
shock should form when supersonic flow from the conjugate 
hemispheres collide at the equator. The flows collide as early as 
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t = 0.25 hour, the flows from the northern hemisphere can be 
seen from Figure 2, and the corresponding flow from the 
southern hemisphere is the mirror image of this 
flow with respect to the equator. When the flow begins to 
overlap, the shock should form through the ion-ion instability. 
The conditions for such an instability in a colliding situation are 
given 

Uv,, SMC,, T,>3T, (8) 

where V ti is the ion thermal velocity, v b is the ion beam 
velocity, C, is the ion-acoustic speed and M is the mach 
number, which could be as large as 4 [Forslund and Shank , 
1970; Montgomery and Joyce y 1969]. However, it must be 
mentioned here that high critical Mach number M = 4 is 
determined by the nonlinear evolution of the electron dynamics, 
including trapping and heating of electrons [e.g., Singh, 1988]. 
For isothermal electrons, as it is assumed in the semikinetic 
model, M= 1.6 [Tidman and Krall , 1971]. It is worth pointing 
out that for large beam velocities, oblique ion waves propagating 
at an angle from the magnetic field are likely to be excited. 
However, the role of such highly oblique waves, which are likely 
to occur f or highly supersonic beams, in momentum exchange 
between interpenetrating beams and shock formation is not well 
understood. 

We examine here the likelihood of the instability occurring 
from the flow parameters given by the semikinetic model. First, 
we do this exercise for t = 30 min, when the flow has crossed 
the equator. Figure 6 shows the average flow velocity V h , the 
temperature ratio T n /T t , and the ion-acoustic speed C, as a 
function of geomagnetic latitude for the flow at / = 30 min., 
shown in Figure 2. Note that the temperature ratio is plotted 
after multiplying it by 10, so that all the line plots in Figure 6 
can utilize the same vertical scale. C, is calculated from 
C t = [k(T'+3T n )lm] V2 ' The critical temperature ratio 
T n / T § = 0.33 for the instability is shown by the segment of the 
thick horizontal line in Figure 6. It is seen that the ions have 
sufficiently cooled down to meet the instability condition on the 
ion temperature over an extended equatorial region (|A| £ 20*) . 
The flow coming from the opposite hemisphere shows a similar 
feature. The ionacoustic speed in the equatorial region is about 
8 km / s ~ l . It is seen that over the latitudinal region (|A| £ 20 s ) , 
the ion beam velocity is about v b s 2 C, . In the semikinetic and 
hydrodynamic models discussed here the electrons are assumed 
to obey the Boltzmann law. Therefore ion beams with such 
velocities are too fast to excite the ion-ion instability and thereby 
to form shocks in the model. Furthermore, it is important to 
point out that the processes which lead to shock formation, 
including the ion-ion instability, are microprocesses, which are 
suppressed in the large-scale models [Singh and Chan , 1993]. If 
electrons dynamics were rigorously included in the model and 
the associated microprocesses properly resolved, it is likely that 
the ion-ion interaction would have occurred forming shocks. 

Figure 7 shows the drift velocity v bl C,, and T„/T, at / = 1 
hour for the flow originating from the northern hemisphere. It is 
seen that as the ion beam penetrates into the opposite 
hemisphere (A<0), it gets progressively warmer and the 
temperature condition 7^/7, £0.3 is not met beyond 
|A| = 10°. Thus ion instability and shock formation are not 
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expected- This indicates that the shock formation in the 
hydrodynamic model (Figures 4 and 5) is an artifact of the 
model. 

In view of the above discussion in connection with Figures 6 
and 7 , it emerges that on the basis of temperature condition 
alone, it can be argued that if shocks form, they should be during 
the early stage when the ion beams begin to interpenetrate in the 
equatorial region. During later stages, when the beams penetrate 
into the opposite hemispheres, the shock formation is not likely 
unless somehow electrons are heated, enhancing the temperature 
ratio T n !T t . However, as mentioned earlier the shock formation 
in the equatorial region requires a rigorous treatment of electron 
dynamics. In view of the simplified treatment of electrons in the 
models described here, and relatively coarse spatial and 
temporal resolutions afforded by them the issue of equatorial 
shock formation cannot be settled in this paper. 

3.3. Subsonic Flow 

After the initial stage of supersonic flows from the conjugate 
ionospheres, the flows become generally subsonic. This is 
predicted from both the models. Figure 8 shows the status of 
the flow at t = 4 hours, from both the hydrodynamic and 
semikinetic models. As before, there are three curves from the 
hydrodynamic model which are compared against the curve from 
the semikinetic model. Figure 8 b shows that the flow velocities 
obtained for different values of rj from the hydrodynamic model 
agree with the flow velocity given by the kinetic model. The 
nunrimum flow velocity of about 5 km / s~ l seen near the 
boundary X = X 0 is subsonic with respect to the ion-acoustic 
speed C, = 10 km / s~ l . We note that the average flow velocity 
peaks slightly above the boundary in both hydrodynamic and 
kinetic models. The peaking is a consequence of the boundary 
conditions at X = X 0 and the acceleration of ions by the pressure 
and electric field distributions in close vicinity of the northern 
boundary of the flux tube. 

Figures 8 a, 8 c, and 8 d show that for the subsonic flow the 
density and temperature structures critically depend on the heat 
flow factor 77 . For rj = 0.3, the hydrodynamic model yields 
results in good agreement with those from the kinetic model. 
When r\ becomes too small (77 = 0.05), the structures in the 
density and temperature profiles markedly differ from the 
kinetic model; the density structure shows an extended density 
cavity in the equatorial region, where parallel temperature is 
relatively high [Singh, 1991]. Furthermore, for low values of rj 
there is density enhancement and correspondingly a low parallel 
temperature in the southern hemisphere. When heat flow factor 
is sufficiently large (7>0.15), such structures in n(X) and 
T^(X) are washed away. In a recent paper. Ho et al. [1993] 
made a similar observation regarding evolution of density 
perturbations in the polar wind by comparing semikinetic and 
hydrodynamic models for an open flux tube. 

The comparison between the hydrodynamic and kinetic 
results at / - 12 hours is shown in Figures 9a to 9d. The density 
and temperature structures at this stage are qualitatively similar 
to that at t = 4 hours as shown in Figures 8 a to 8 d. However it 
is seen that at t - 12 hours, the density and temperature profiles 
even for 77 = 0.05 have begun to compare well with that for 
77 = 0.3, for which the density distribution agrees well with that 
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given by the semildnetic treatment The discrepancy between the 
densities predicted by the kinetic treatment and the 
hydrodynamic one for rj — 0.3 is bounded by 1 5%, for most part 
of the flux tube, except near the southern boundary x - -X 0 . 

The runs for the comparison between the hydrodynamic and 
kinetic models were carried on until t - 48 hours. For time 
t > 12 hours, it was found that the hydrodynamic model 
systematically yields densities higher than that given by the 
semildnetic model. Figure 10 shows the comparison between 
the equatorial densities as obtained from the two models. This 
figure shows the te mp oral evolution of the equatorial densities 
found from the kinetic (solid line) and the hydrodynamic 
(dashed lin g curves) models. For the latter model, the densities 
are plotted for different values of the flow velocity V 0 at the 
boundaries X = ±X 0 . We remind ourselves that the results from 
the hydrodynamic model shown in Figures 3, 4, 5, and 8 are for 
a flow velocity V 0 = (kT 0 1 2jvn) vn - 0.39F, . We notice from 
Figure 10 that for this boundary value of the flow velocity, the 
kinetic and hydrodynamic curves are remarkably close for f = 12 
hours. This implies that this boundary value of the flow velocity 
closely corresponds to the input flux determined by a half- 
Maxweilian distribution function, which is imposed as a 
boundary condition in the kinetic model. For t > 12 hours, the 
boundary value of V a = 0.39F, yields an overrefilling compared 
to the kinetic model. This simply implies that the net influx of 
ion into the flux tube at the ionospheric boundaries steadily 
decreases in the kinetic model, primarily due to the ions flowing 
out of the flux tube. On the other hand, in the hydrodynamic 
model, the influx is primarily determined by the imposed flow 
velocity, and it remains constant This is demonstrated by 
comparing the temporal evolution of the total plasma content in 
the flux tube as seen from the two models. 

Figure 1 1 shows the total content as a function of time. As in 
Figure 10, for the hydrodynamic model the curves are for 
different values of the imposed velocity at the boundary. It is 
seen that for V Q - 0.39F,, the hydrodynamic model yields nearly 
the gflnv total plasma content as the kinetic model with nearly 
the same rate of increase in it for / = 12 hours. At later times, 
the content from the kinetic model shows a tendency toward 
saturation because the rate of increase in the content 
continuously decreases. Even in the hydrodynamic model there 
is a tendency toward the decreasing rate, but the decrease is 
much slower. This difference in the influx of the ions from the 
two models has a simple explanation. In the kinetic model, some 
of the ions have the liberty to exit the flux tube as they are 
scattered by Coulomb collisions, or as they simply flow out On 
the other hand, in the hydrodynamic model the plasma entering 
the flux tube can leave the system only through the opposite 
boundary, where the flow velocity becomes exceedingly small 
after the shock phase ( t > 2 hours). This implies that in the 
hydrodynamic model there is no provision for the plasma to 
leave the system. The slight tendency toward the saturation in 
the hydrodynamic model is due to the changing plasma condition 
near the boundary where the flow originates. As the plasma 
density near this boundary increases, the influx into the flux 
tends to decrease. 

Figures 10 and 1 1 also show the equatorial density and the 
total plasma content from the hydrodynamic model for V 0 = 0.1F, 
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and V 0 =0. It is seen that even or v o = 0, the equatorial density 
and the total content are increasing with time and, in about 48 
hours, they tend to approach the corresponding results from the 
kinetic modeL It may 9 ound strange how refilling can occur with 
a boundary condition of V 0 — 0 ! The refilling of a flux tube with 
zero flow velocity as boundary conditions in a hydrodynamic 
model was previously described by Singh et al. y [1986]. 

The ion flux developing near the boundaries from the kinetic 
model and the hydrodynamic model for V 0 = 0.39F, during the 
early firngs (f < 12 hours) is L8 x 10® ions cm~ 2 s~ l . In the kinetic 
model, this flux continuously decreases because some of the ions 
entering the flux tube eventually leave. However, in the 
hydrodynamic model the plasma entering the flux tube remains 
in it, causing the ovorefilling for t > 12 horns when 
v 0 ** 0.39F, . When „ =0, the hydrodynamic model yields a flux 
of £10* ions , and it continuously decreases as the 

forces (determined by density and temperature gradients) on the 
ions accelerating them into the flux tube from the boundary cells 
diminish es with the refilling. It is worth mentioning that the 
comparison carried out above is based on the simplified 
boundary conditions in the kinetic and fluid treatments and a 
heuristic treatment of the heat flow. A comparison of the plasma 
treatments without these simplifications will be worthwhile. 

A comparison of plasma distributions in the flux tube at a 
relatively late time ( / = 48 hours), as obtained from the two 
models, is shown in Figures 12a to 12cL For the hydrodynamic 
model, 7 = 0 . 3 , and the distributions are given for three values 
of the boundary velocity; V Q = 0.39F r 0.1F,, and 0. Density 
profiles in Figure 12a show the ovenefilling for V 9 = 0.39F ; , 
but when V 9 is reduced below 0. 1 V t , , the density profiles from 
the two models disagree near the boundary k = k ot but away 
from it the agreement considerably improves. The disagreement 
near the boundary is also reflected in the velocity profiles in 
Figure 12b. Despite the above disagreement in the density and 
velocity profiles near the boundary, the temperature structures 
obtained from the two models are nearly identical. Temperature 
is nearly isotropic it rises to about 0.4 eV in the 

equatorial region from the boundary value of 0.3 eV. In the late 
stage of the refilling, the similarity between the temperanne 
profiles, despite the differences in the density and velocity 
profiles from the two models, can be understood by examining 
the temperature equations (3) and (4) and the flow properties. 
Since the flow velocity is small, the velocity profile has little 
effect on the temperature profiles. The maximum flow velocity 
in a localized region near the boundary is 2.5 kms~ l , compared 
to the thermal velocity of 5.5 kms~ [ and ion-acoustic speed of 
10 kms~ l . The density distribution affects the temperature 
distribution through the heat flow terms in equations (3) and (4). 
In the late of the flow when the gradients have smoothed 
out and the densities are relatively large, the density distribution 
also has insignificant effects on the heat flow. 

4. Conclusion and discussion 

We have carried out a comparison between semikinetic and 
hydrodynamic models for plasma flow along closed magnetic 
field lines. The comparison has direct relevance to the problem 
of plasmasphehc refilling. It is found that the comparison does 
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not depend only on the plasma physics afforded by the models, 
but it also strongly depends on the boundary condition on the 
flow velocity. In a kinetic model, an appropriate boundary 
condition is to prescribe the velocity distributions of the 
inflowing ions to be haif-Maxwellian for V £ 0 at X = X 0 and 
for at X - -X 0 . In the hydrodynamic model this 

boundary condition corresponds to the drift velocity 
V 0 =(kT 0 f2xm) l/ 2 . A comparison of results from the two 
models with such boundary conditions revealed the following 
important features of the flows. 

1. When supersonic flows develop in response to a sudden 
depletion in a flux tube, the hydrodynamic and kinetic models 
yield distribution of density, flow velocity, and temperatures in 
generally good agreement The temperature distributions in the 
region of supersonic flows are found to be remarkably similar, 
showing small effect of the heat flow. It is worth pointing out 
that Deman and Shunk [1991] compared the behavior of a 
highly supersonic plasma flow from a hydrodynamic model 
based on a more complete (16-moment) set of equations with 
that from a semikmetic model, demonstrating a good agreement 
We have demonstrated here that, for a highly supersonic flow, 
even a much simpler set of hydrodynamic equations are 
adequate. It is physically explained by the fact that the transport 
of heat in a supersonic flow is dominated by the large drift 
velocity and not by the heat flow process. Mathematically 
speaking, it implies that the heat flow terms are negligibly small 
compared to the convective terms in the temperature equations. 

2. Both models show reflection of the supersonic flow when 
it penetrates deep into the opposite hemisphere. Since even a 
two-stream hydrodynamic model cannot handle the 
counterstreaming for a given flow, the reflection automatically 
leads to a shock formation [Rasmussen and Schunk, 1 988; Singh, 
1991]. The shock first moves upward toward the equator and 
then downward to the ionospheric boundary. An examination of 
the plasma conditions for shock formation shows that the shock 
seen in the hydrodynamic model is an artifact of the model; the 
ion beams are found to be too warm to excite the ion-ion 
instability which can subsequently produce a shock. The 
semikinetic model shows the development of counterstreaming 
for the flow; the counterstreaming advances to the equator and 
downward to ionospheric boundary. It turns out that the transit 
time of the shock all the way to the ionospheric boundary and 
the time for the counterstreaming to spread to this boundary are 
nearly the same, about 2 hours. In a previous paper, Singh 
[1991] reported the shock transit time to be about 4 hours, which 
is in error due to a normalization factor of 2. In view of the short 
transit time of the shock, the shock formation does not 
significantly affect the refilling, as evidenced by the comparison 
of the flows from the hydrodynamic and kinetic models for later 
times. 

Lack of shock formation in the equatorial region, when the 
ion beams begin to interpenetrate [Banks et al., 1971] is 
uncertain in view of the spatial and temporal resolutions 
afforded by a large-scale model and the simplicity in handling 
the electron dynamics by the Boltzmann Law. 

3. After about 2 hours, the flow in each hemisphere becomes 
subsonic with respect to the ion-acoustic speed. This is seen 
from both models. 
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4. A comparison of the total plasma contents and the 
equatorial densities from the two models indicates a good 
agreement up to about t = \2 hours, after which the 
hydrodynamic model indicates overrefillmg of the flux tube. The 
overrefilling is traced to the inability of our hydrodynamic model 
to control the net p lasma inflow by the returning particles. The 
inflow is determined by the imposed boundary conditions and 
the outflow of plasma is exceedingly small. On the other hand, 
in the kinetic model the influx gradually decreases due to ions 
returning from the flux tube, showing a tendency toward 
saturation in the refilling in about 2 days. It is worth pointing 
out that it will be useful to perform a study comparing the 
models based on the kinetic and hydrodynamic treatments by 
relaxing some of the simplifications in terms of boundary 
conditions and in handling the heat flow in the latter treatment 
The boundary conditions can be relaxed by including the 
ionospheric plasma generation processes at low altitudes [Gutter 
and Gombosi, 1990]. 

When the boundary flow velocity in the hydrodynamic model 
reduces below V 0 =(kT 0 I2xm ) U2 , there is an initial 
underfilling, but eventually, the refilling from this model catches 
up to that given by the semikinetic model. For example, when 
V a 0, the degree of refilling from the two models, in terms of 
both the equa tori al density and the total plasma content in the 
flux tube, becomes approximately the same in about 2 days. 

In some previous studies [ Singh et al 1986; Rasmussen and 
Schunk, 1988; Singh, 1991], a boundary condition of zero flow 
velocity was used. It may appear strange that a refilling occurs 
with thia boundary condition on the flow velocity. The issue is 
briefly revisited here. 

Fran the comparison of the plasma contents and the 
equatorial densities given by the models, it is concluded that 
after about 12 hours, the choice of boundary condition in the 
hydrodynamic model is quite uncertain. In view of this 
uncertainty, the choice of zero-velocity boundary could be useful 
during the late stage of the refilling, it yields undemefilling only 
near the boundaries, where the density and average flow velocity 
show a discontinuity in the flow, otherwise, over the rest of the 
flux tube the density and flow velocity are in quite good 
agreement with those given by the kinetic model. 

The hydrodynamic model described here is a two-stream 
model and includes equation for the parallel and perpendicular 
temperatures. Single-stream hydrodynamic models [Singh et al., 
1986; Gutter and Gombosi , 1990] suffer from the shortcoming 
that they generate shocks at the equator whether the plasma 
conditions allow them or not The single- and two-stream 
models with assumed temperature isotropy suffer from the 
shortcoming that the shock transit time is fairly long, and a 
major part of the refilling occurs through supersonic flows from 
the ionospheres ( Singh et al., 1986; Rasmussen and Schunk , 
1988; Singh, 1991]. This ir in contrast to the two-stream model 
with a self-consistent treatment of the temperature anisotropy, 
this model yields evolution from supersonic to subsonic flows at 
the same timescale as the kinetic model. In this sense the 
heuristic treatment of the heat flow described in this paper 
appears to be adequate. This treatment also appears to be 
adequate even in the subsonic stage as long as the flow velocity 
near the boundary is relatively large near the thermal speed, for 
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example, for f < 12 hours in our present calculations. However, 
when the flow velocity becomes sufficiently low so that a large 
fraction of injected ions in the kinetic model begin to return 
from the immediate vicinity of the boundary, the boundary 
conditions for the hydrodynamic model diverges from that of the 
kinetic treatment because this model does not allow for a return 
flux for a given stream. Can this situation be improved by a 
more rigorous treatment of the heat flow and/or by properly 
including the ionospheric plasma supply [Gutter and Gombosi, 
1990]? In order to answer this question, it will be useful to 
compare models based on (1) the heuristic heat flow 
treatment, (2) a more sophisticated treatment of the heat flow 
using 16-moment set of transport equations, and (3) the 
semikinetic treatment, and all models properly including the 
analogous ionospheric boundary conditions. If the model (1) 
compares well with the latter ones, the computational effort in 
using transport equations for modeling space plasma will be 
considerably reduced. 
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Figure 1. Geometry of a closed magnetic flux tube. The 
latitudinal angle X and the geocentric distance r are shown. The 
ionospheric boundaries are at s = 0(X = X 0 ) t and 


Figure 2. Temporal evolution of the flow originating from the 
northern hemisphere; phase-space (s-V n ) plots are shown. The 
plot at / = .003 hour nearly shows the initial plasma in the 
flux tube. 

Figure 3. Comparison of flows from semikinetic and 
hydrodynamic models at f = 30 minutes. For the latter model 
flows are given for 7 = 0 , 0.05 and 0.3: (a) density, (b) flow 
velocity, (c) parallel temperature, and (d) perpendicular 
temperature distributions. For most latitudes the curves are so 
close together that it is difficult to distinguish them. 
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Figure 4 . Same as Figure 3, but at t = 1 hour. The 
hydrodynamic curves are distinguished by the presence of a 
shock which is nmni frated by sudden jumps in density, flow 
velocity and parallel temperature near X s -25. Shocks are 
indicated by the arrows. 

Figure 5. Propagation of the shock is shown through the jump 
in the density profiles at (a) / = 1 hour, (b) t =1.5 hours, and 
(c) t = 2 hours. For the purpose of comparison the curve from 
the kinetic model and three curves for ij = 0, 0.05, and 0.3 from 
the hydrodynamic model are shown. Shocks are indicated by 
arrows. 

Figure 6 . Distribution of flow parameters from the kinetic 
model. The average flow velocity V b , ion-acoustic speed C s and 
the temperature ratio 7 J, / 7 ^ for ions with > 0 are shown for 
the purpose of instability analysis at t = 30 minutes. The 
corresponding curves for F [,<0 can be deduced from the 
symmetry considerations. 

Figure 7. Same as Figure 6 , but at r = l hour. 

Figure 8 . Same as Figure 4, but at t = 4 hours. 

Figure 9. Same as Figure 4, but at t = 12 hours. 

Figure 10. Temporal evaluation of the equatorial density from 
the kinetic model and from the hydrodynamic model for three 
values of V ot 0.39 V t ,0.1 V Jt and 0V r 

Figure 11. Temporal evolution of the total plasma content in 
the flux tube for the kinetic and hydrodynamic models. 

Figure IX Distribution of (a) density, (b) flow velocity, (c) 
parallel temperatures and (d) perpendicular temperature at t = 
48 hours. All vertical scales arc linear in this Figure. There are 
three curves from the hydrodynamic model for F o = 0, 0.1 V t , 
and 0.39 V f . 


Figure 1 . Geometry of a closed magnetic flux tube. The latitudinal angle X and the geocentric distance r are 
shown. The ionospheric boundaries are at j = 0 (A = X Q ) , and jr = s WtMX (X = X 0 ). 

Figure X Temporal evolution of the flow originating from the northern hemisphere; phase-space (s - V^) plots 
are shown. The plot at / = 003 hour nearly shows the initial plasma in the flux tube. 

Figure 3. Comparison of flows from semikinetic and hydrodynamic models at t = 30 minutes. For the latter 
model flows are given for 7 = 0, 0.05 and 0.3: (a) density, (b) flow velocity, (c) parallel temperature, and (d) 
perpendicular temperature distributions. For most latitudes the curves are so close together that it is difficult to 
distinguish them. 

Figure 4 . Same as Figure 3, but at t = 1 hour. The hydrodynamic curves are distinguished by the presence of 
a shock which is manifested by sudden jumps in density, flow velocity and parallel temperature near X = -25. 
Shocks are indicated by the arrows. 
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Figure 5. Propagation of the shock is shown through the jump in the density profiles at (a) / — 1 hour, (b) 
r =1.5 hours, and (c) f - 2 hours. For the purpose of comparison the curve from the kinetic model and three 
curves for rj = 0, 0.05, and 0.3 from the hydrodynamic model are shown. Shocks are indicated by arrows. 

Figure 6. Distribution of flow parameters from the kinetic model. The average flow velocity V b , ion-acoustic 
speed C, and the temperature ratio T n / T t for ions with v n > 0 are shown for the purpose of instability analysis 
at f = 30 minutes. The corresponding curves for ^, < 0 can be deduced from the symmetry considerations. 

Figure 7. Same as Figure 6, but at r = 1 hour. 

Figure 8. Same as Figure 4, but at t = 4 hours. 

Figure 9. Same as Figure 4, but at / = 1 2 hours. 

Figure 10. Temporal evaluation of the equatorial density from the kinetic model and from the hydrodynamic 
model for three values of V Qt 0.39 V n 0.1 V /9 and OF,. 

Figure 11. Temporal evolution of the total plasma content in the flux tube for the kinetic and hydrodynamic 
models. 

Figure 12. Distribution of (a) density, (b) flow velocity, (c) parallel temperatures and (d) perpendicular 
temperature at / - 48 hours. All vertical scales are linear in this Figure. There are three curves from the 
hydrodynamic model for ^ = 0,0.1^, and 0.39^. 
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AbStraC L piasmaspheric refilling have revealed that during the early stage of the refilling 

i^arge muu f ^ . u/wvpwr the instability of such ion beams and its effect on 

in'S?we weTeat ion kineiically and electrons are assumed to obey the Boltzmann law In the otter 

SHrr^ 

SsssS-ssr 

enhanced plasma fluctua f t0 excite ^ instability. The former instability heats the 

double layers greatly affecting the state of the plasma inside the central region of the flux tube. 

1. Intro ^ t ‘®" tandi prob i em i n S p a ce plasma transport is the coupling of microscale and mesoscale Presses. 

i nf thic nrnhlem is the early stage refilling of the outer piasmaspheric flux tubes after 

S2rS si-ss ^ =2= 

LerhcmSpheric plasma flow. The counterstreaming ion beams should be able to drive instabilities and affect th 

pusnrnfl^^d te^mo mmng pr^ ^ ^ dliri „ g the ,nst ton years or so have been 

primarily^ lSg™,e modeling using hydrodynamic [tennnv « cl. 1984; So.** e, cl. 986; ««» W 
y 1Q o 2 . sineft 1988 1990 1991a; Singh and Chan, 1992] and semikinetic [Wdson et al, 1992, Lin et al, 
% 92 ] ffeatmente The temporal and spatial resolutions afforded by the model are generally too crude to include 
the microscooic effects Driven by the idea that interpenetrating ion beams may generate shocks during refilling 
Lfa «Tl9711 the physios of shock formation has been studied using small-scale stmulauons [Sug/l and 
[Banks ec at i l, , i q«« ^ n gh and Chan 1992]. More recently some anomalous effects m 

EatSe is caused by waves generated by a hot plasma trapped in the equatorial region [Stngh et al., 1982, Singh 
scattering is cau^a oy s Singh 1991b; Singh and Chan, 1992; Lin et al., 1992], 

tostndy the effect of iondnmm driven wrves on rite rilling of a magneric 

(tax tnbe™ri.7arSowing into iUtom its ends, as schematically shown in Figure ia. Stnce we wtsh o revive 

this question in this paper. 
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The basic ideas found in this paper are as follows. When the simulation is performed by treating ions 
kinetically and electrons by simply assuming that they obey the Boltzmann law, the ion-beam driven instabilities 
are not found because the beams are too fast to excite the waves. On the other hand, when electrons are also 
treated kinetically the fast ion beams during the early stage of the evolution of the plasma in the flux tube drive 
electron-ion instability [Fujita et al., 1977; Singh, 1978], This instability contributes to the heating of electrons in 
conjunction with the electron energization in the self-consistent potential distributions resulting from 
counterstreaming plasma expansion. This heating increases the ion-acoustic speed in the plasma and thereby the 
initially fast ion beams eventually become slow enough relative to the ion-acoustic speed to drive the ion-ion 
instability, which is found to be very effective in coupling counterstreaming ion beams and in trapping of plasma in 
the flux tube, especially in its outer regions away from the central minimum magnetic field (Figure lb). 
Comparing the results from the simulations with and without the kinetic treatment of the electrons, we found that 
in the former simulation the filling of the flux tube is enhanced due to the trapping of plasma in the potential 
structures set up by the ion-ion instability. The instability also drives electrostatic shocks and distributes the 
plasma in the flux tube quite differently from that in the simulation with the simplified electron dynamics. 
Furthermore, this instability controls the influx of ions into the central region of the flux tube, and thereby the 
dynamics of the potential structures including the electrostatic shocks in the central region. 

The rest of the paper is planned as follows: In Section 2 we describe the numerical technique. In Section 
3 numerical results are described along with their interpretation. The conclusion of the paper is given in Section 4. 

2. Numerical Model 

We perform a one-dimensional particle-in-cell simulation of plasma flow along an artificial flux tube 
(Figure la). The magnetic field B(X) = B a {\- a exp[-(X-i//2) 2 /o 2 ]) where B 0 is a constant field 

outside the minimum-field region (Figure lb), d is the size of the simulation system and the choice of a and a 
determines the desired field distribution. The cold plasma flows into the flux tube from the two plasma reservoirs 
at the ends of the simulation system at X = 0 and X = d (Figure la). The simulation technique is described in 
Singh and Chan [1992]. Plasma dynamics is simulated using a particle-in-cell (PIC) code. We have performed 
two types of simulations; in one type, called here Run-1, ions are treated kinetically using the PIC code while 
electrons are assumed to obey the Boltzmann law. This assumption along with the assumed condition of a 
quasineutral plasma gives the electric field E = ~(JcT 0 / ety) dt\ / dx where is the ion density, k is the 
Boltzmann constant, and T a is the electron and ion temperatures in the plasma reservoirs. In the other type of 
simulation, called Run-2 here, both electrons and ions are treated kinetically using the PIC code and we solve the 
Poisson equation for the electric field with the boundary conditions <p(x -0) = <p(x = d) = 0, where <p(x) 
denotes the potential distribution. As the particles move in the flux tube, their magnetic moments are assumed to 
be conserved. 

The velocity distribution function of the ions and electrons in the plasma reservoir is assumed to be 
Maxwellian with a temperature T 0 . The reservoirs supply a continuous flux of charged particles into the flux tube 
through the process of plasma expansion. In the simulations reported here, we have used m i Jm e = 400, which 
adequately separates the electron and ion time scales and, at the same time, allows computationally feasible runs. 
Numerical parameters of the simulations are as follows: system size d = Sx \Qp Z d , the magnetic field 

parameters a = 0.9, and a = 150 Z d for which B(x) is plotted in Figure lb, cell size Ax = 20 Zj, and 
time step At = 0. 1 0)~^ eo , where Zj and (O peo are the Debye length and electron-plasma frequency in the 

cold plasma reservoirs, respectively. In the following discussion we have used normalized quantities defined as 
follows: time /=/ (O pio , distance X = X / Zj, velocity V = V / V te and electric potential 0=e0/kT o , 

where V te = ( kT a / )'~ , (Q pio = (m e / <O peo , and (O peo is the electron plasma frequency in the 

plasma reservoirs and the Debye length (Opeo- 
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3. Flux Tube Filling and Ion-Beam Driven Instabilities 

In this section we examine the evolution of the plasma and the potential distributions in the flux tube as 
seen from Run-1 and Run-2. A comparison of results from these runs reveal when and why plasma instabilities are 
excited and how they affect the plasma and field distributions. 

Figure 2 shows the evolution of the ion phase space in X -Vj\\ plane from Run-1. Note that the velocity 

on the vertical axis is measured in units of the ion thermal velocity V ti = ( kT a / Figure 2a for 7 = 200 

shows the expanding ion beams from the two plasma sources shown in Figure la. By the time 7 = 10 3 (Figure 

2b) the counterstreaming of ion beams is set up all along the flux tube; the counterstreaming continues without a 
significant interaction between the beams as seen from Figures 2c and 2d for t - 2000 and 4000, 
respectively. The main feature of the spatial evolution of the ion beams is that they progressively slow down as 
they approach the opposite end of the flux tube. From the plot at 7 = 4000, we also note that a few particles in 
the ion beams approaching the opposite ends slow down upto almost zero velocity, the number of such ions 
increases with time as seen from the plots at / = 8000. This tends to slowly fill in the velocity space between the 
two counterstreaming ion beams in the off-central region. We find that such scattering of ions is not caused by 
anbipolar electric fields, which are given by the ion density gradients. 

The evolution of the ion density profile in the flux tube is shown in Figure 3. After the initial stage 
(/ <1000), the density profiles are generally smooth having relatively large gradients in the control region 

(1500 < \X\ < 2500), and in the off-central region (100 < X < 1 500 and 3500 < X < 4900) the gradients 

are weak giving weak electric fields. The scattered particles in Figures 2d and 2e occur in the latter regions. In 
view of the weak anbipolar electric fields, incapable of slowing down or reflecting the beams, we suggest that the 
scattering occurs via some anomalous effects involving fluctuations in the plasma. Later we discuss this issue 
further after we have considered the flow in Run-2. 

The evolution of ion phase space in X-Vj\\ plane for Run-2 is shown in Figures 4a to 4f. Note that in 
this figure the velocity is normalized with respect to , in contrast to in Figure 2 for Run- 1 . A comparison 

of the phase-space plots in Figure 2b and Figure 4b, which are for the same time 7 = 10 3 , shows a significant 

difference in the topology of the distribution of the average flow velocity of ions in the two runs; in Run-1 with the 
Rnftymann electrons, the ions have relatively large velocities near X = X mi j creating a "bulge" in the plots. 
The bulge is not seen in Run-2. The bulge in Figure 2 is simply a manifestation of the magnetic field distribution 
B( X) (Figure lb), which controls the density distribution and hence the potential distribution according to the 
Boltzmann law. The relatively sharp density gradients in the density distribution on either side of the midpoint 
yield relatively large electric fields, which accelerate ions coming into the central region and then decelerates them 
while leaving. On the other hand, in Run-2 the Poisson equation including a self-consistent treatment of the space 
charges does not yield a potential distribution controlled by X') . The potential distributions for Run-2 are 

shown in Figures 4g to \t . 

In Run-1 we saw only a weak scattering of the beam ions at late times (l > 4000) into the velocity 
region between the two counterstreaming beams outside the central region of the simulation region. In contrast, in 
Run-2 we find that wave-particle interactions plays a significant role in coupling and mixing the counterstreaming 
ion beams. This can be seen from the phase-space plots for 7 > 1400 in Figures 4c to 4f. The plot at / = 1400 
shows that ion-ion instability has occurred outside the central region and it is sufficiently strong to locally couple 
the ion beams. The instability is more clearly manifested in the potential plots given in Figures 4g to 4 L . The 
corresponding electron phase space plots are shown in Figures 4m to 4r. Both electron and ion phase space show 
formation of vortices, which evolve from the ion-ion instability driven by the counterstreaming ions. We now 
examine in detail why this instability is so prominent in Run-2, but not in Run-1. 

Ion-Ion Instability: The conditions for the ion-ion instability in a plasma having counterstreaming ion beams 
depends on a number of parameters of the beam-plasma system: these parameters are the relative beam velocity 
Vf, r between the two beams, the ratios of the beam temperatures T b t and T b2 to the electron temperature T e , 
and the relative densities of the ion beams fi b \ / W>"ied and Wong, 1986]. The relative velocity and the 
temperature ratios T bX / T e and T b2 I T e play crucial roles in determining whether the instability occurs or not. 
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For two symmetric counterstreaming ion beams, i.e., V br = 2V b , where ±V b are the two beam velocities, and 
^b\ ~ ?b 2 ~ Ti , the instability occurs when [McKee, 1970] 

l3V ti <V b <C s and T e > 32J (1) 

where C s is the ion-acoustic speed given by C s = (k( T e + 3 T^) / rr^ ) 1/2 and 

For non-symmetric beams it is not possible to give analytical conditions as above for the symmetric ones. 
However, some qualitative assessments can be made; if one of the beams is warmer than the other, say ^1 > ^2’ 
instability occurs as long as the velocity condition in equation (1) is met and T e / T b j > 3 [Baker, 1973; Gresillon 
and Doveil , 1975; Singh, 1978]. However, Fried and Wong [1966] have shown that when one of the beams is 
sufficiently cold depending on the relative beam density, the instability can occur even if the warm beam satisfies 
the condition T b \f T e = 1 . 

Now we examine if the instability conditions described above are met in our simulations. Figure 5 shows 
examples of ion velocity distribution function from Run-1; the distributions shown here are from the off-central 
region 4000 < X < 4500, where relatively slow ion beams are more likely to meet the condition for the ion-ion 
instability. Note that in this figure the horizontal scale is the ion velocity normalized with respect to the thermal 

velocity V ti = (k T Q lm$^. The distributions shown in Figure 5 are for (a) /=4000, (b) /= 6000, and (c) 
t = 8000; we noticed earlier from Figure 2 that at such late times the beam ions weakly scatter to fill in the 
velocity space between the two ion beams in the off-central region. 

The instability conditions discussed above are for a Maxwellian velocity distribution, but we note from 
Figure 5 that the distributions are not Maxwellian. However, effective beam velocities and temperatures as first 
and second order moments of the individual ion beams can be calculated, and they are given in Table 1. Note that 
the velocities given in this table are normalized with respect to V (i . The ion beams for V\\ < 0 and V\\ > 0 are 
called beam-1 and beam-2, respectively. The relative beam velocity V br — and the ion-acoustic speed 

are also tabulated. We find from Table 1 that the beam temperatures do satisfy the condition for the instability 
given by equation (1). However, the relative beam velocity is too large to excite the instability, that is, V br > 2 C,; 
thus the ion-ion instability, as predicted from the linear instability analysis for Maxwellian ion beams, should not 
occur. However, as noted earlier, the beams are not Maxwellian; the fall of the distribution functions at velocities 
below their respective peak velocities is quite sharp and the peak separation in velocity progressively decreases 
from AV =3.1 at / =4000 to AV =2.6 at t - 8000. Taking AV as V br we find that the 
instability condition on velocity is marginally satisfied at t = 8000. The approach towards the marginal 
instability causes an enhanced fluctuation in the plasma (e.g., see Ichimaru [1973]) and the consequent scattering 
of ions in velocity space as mentioned earlier. Figure 6 shows evolution of the density fluctuations at 
X = 1 500, 1000, and 500; we note from this figure that with the decreasing distance from the end of the flux 
tube, the amplitude of the density fluctuation increases and therefore the fluctuations in the electric field also 
increases. The fluctuation in the electric field scatters the ions, populating the velocity region between the ion 
beams. The scattered ions appear as a "tail" to the velocity distribution function near V = 0 as seen from Figure 
5c. In the central region of the flux tube, where the ion beams are quite fast and far from the marginal instability 
condition, the fluctuations in the density and the electric fields are weak and scattering of ions is not seen (see 
Figure 2e). 


Table 1. Ion beam parameters from Run-1. 


t 

n b L 


Vb\ 


7 br 

%X 

%2 

T 

H 

4000 

0.55 

0.45 

-1.60 

1.85 

3.45 

0.16 

0.18 

l 

1.22 

6000 

0.52 

0.48 

-1.58 

1.60 

3.18 

0.2 

0.23 

l 

1.30 

8000 

0.51 

0.49 

-1.51 

1.51 

3.00 

0.21 

0.27 

l 

1.31 
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The evolution of the ion and electron velocity distribution functions from Run-2 are shown in Figure 7; 
Figures 7a, 7b and 7c show the ion distribution functions at t = 1000, 1200 and 1400, respectively, and the 
corresponding electron distribution functions are shown in Figures 7d, 7e and 7f. The distributions shown in these 
figures are for the plasma in the region 4000 < X < 4500 and for times just before the onset of the ion-ion 
instability leading to the formation of vortices in Figure 4. Table 2 shows the ion beam and plasma parameters 
relevant for the linear instability analysis. The ion beams denoted by subscripts T and "2" are for V\\ < 0 and 
K > 0, respectively. All velocities given in this table are normalized with respect to Vie- Comparing 
v br and C s in the table, we find that by the time / = 1400, the velocity condition on the instability is met, 
i.e., V br <2 C s . We note that the temperature condition for the i-i instability is also satisfied because 

T bi / T e = 0. 13 and T b2 / T e = 0.3. Thus, the instability occurs as manifested by the ion and electron phase- 
space plots, and the potential distribution in Figure 4 for times t > 1400 . 

Table 2. Ion beam and plasma parameters from Run-2. 


t 

m 

1^ 

wm 


EH 

mm 

% 2 

T g 

mm 

1000 

0.93 

0.07 

-0.1 

0.22 

KB 

0.32 

0.58 

2.6 

0.1 

1200 

0.87 

0.13 

-0.096 

0.19 

BUB 

0.37 

0.74 

2.8 

0.1 

1400 

0.82 

0.18 

-0.086 

0.16 

0.21 

0.42 

1.00 

3.1 

0.12 


Electron-Ion Instability: It is noteworthy that in Run-2 electron heating is the key for fulfilling the instability 
condition by increasing the ion-acoustic speed; Table 2 shows that the electrons are heated to an effective 
temperature of 3.1 T a by t = 1400. How are the electrons heated? We discuss here the mechanism for the 
electron heating. 

In an ion beam-plasma system, e-i instability is yet another mechanism for exciting waves, which heats 
the electrons [McKee, 1970], The e-i instability occurs when the ion beams are relatively fast [Fujita et at., 1977; 
Singh, 1978], This instability is kinetic and it occurs as the slow (negative energy) wave of an ion beam undergoes 
Landau damping by the electrons [Hasegawa, 1975; Singh 1978]. Such an instability creates density fluctuation at 
frequencies less than the ion plasma frequency. In order to demonstrate that this instability does occur in Run-2, 
we show in Figure 8a and 8b the evolution of the plasma and the potential by plotting n e (t), «■(/) and <p(t) 
at X = 2500 . Figure 8a shows the evolution of the electron and ion densities. For t < 200, the central region 
is nearly unpopulated with a plasma because the expanding plasmas from the sources at the end of the flux tube 
have not yet reached there. From t ~ 200 to / ~ 400, the central region is electron rich, and when t > 400 
the central region acquires a quasineutral plasma with n e = /J, . The plasma density builds up to 

~10 -2 n o by f = 1000. 

The electric potential shown in Figure 8b is expected to slowly evolve with the plasma build up in 
response to the counterstreaming plasma expansion into the flux tube. This slow evolution is shown by a dotted- 
line curve, which is the potential averaged over a time interval At = 100. The fluctuations superimposed on this 
curve, as shown by the solid-line curve, are relevant to our discussion here. The frequency spectrum of the 
potential variation in Figure 8b is shown in Figure 8c. We notice that there is a broad peak in the frequency range 
0.01 < / / f pio <0.1, where f pio is the ion plasma frequency in the source plasmas. For the time period of 

simulation shown in Figure 8, the local plasma density at X = 2500 is n <10 ~ n n , and hence the local ion 
plasma frequency f pi <0.1 f pio . Thus the low frequency oscillations are bounded by this local plasma 
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frequency. These low frequency fluctuations which have amplitudes upto 4 kT a l e , heat the electrons to a 
temperature of upto 3 T a as shown in Table 2. 

Another mechanism which contributes to the electron energization deals with the spatial distribution of 
the slowly time-varying potential as shown in Figure 4g. When the electrons expanding from the plasma sources 
cross the midpoint of the simulation region, they see an accelerating potential. A signature of such acceleration is 
the elongated tail of the electron velocity distribution function for > 0 at t = 1000 shown in Figure 7d. 

Consequences of the Ion-Ion Instability: The plasma flow and the potential distribution in the flux tube are 
significantly affected by the ion-ion instability as seen from Figures 4a to 4r. The instability occurs in the off- 
central region and in the central region the ion beams continue to counterstream (Figures 4d to 4f). At an early 
stage, the instability occurs at relatively small scale lengths as seen from the potential distribution in Figure 4i. 
The subsequent nonlinear evolution of the instability forms vortices, which grow in size as clearly seen from the 
potential distributions for t ^ 1 600 in Figures 4j, 4k and 4 ( . The amplitudes of the potential associated with 
the vortices are as large as 5 kT a / e, which approximately corresponds to the ion beam energy at t = 1400 
when the instability sets in, that is, 

where A<f> is the amplitude of a vortex and is the beam velocity of the faster ion beam. The above 
relationship between V b 2 and A(f> suggests the possibility of trapping the ions in the vortices. This trapping is 
evident from Figure 4. The cells with relatively high positive potentials have ions almost completely mixed in 
velocity space and their velocity spread is reduced in contrast to the relatively large velocity spread in the 
neighboring negative potential cells. As a matter of fact in the latter cells ions have a M hole H in their velocity 
distribution function, which is a signature of an ion hole. The electrons 1 behavior is directly opposite to that of the 
ions; electrons have large velocity spreads in the positive potential cells and a relatively small spread in the 
negative potential cells (see Figures 4p, 4q and 4r). This alternate trapping of electrons and ions affects the flow of 
plasma and density distribution in the flux tube. Before we discuss the density distribution and the filling of the 
flux tube, we point out that another consequence of the ion-ion instability is the formation of electrostatic shock- 
like structures. 

Electrostatic Shocks: The electrostatic shock-like structures evolve from the ion-ion instability; the vortices 
having positive potential cells begin to grow in size with increasing time and their edges steepen into shock pairs. 
The potential distribution plots for t > 1400 show examples of the formation of some shock pairs as indicated by 
arrows in Figures 4j, 4k and 4 1 . Inside a shock pair the potential is positive; ions are mixed and electrons are 
trapped showing their acceleration above the background electrons. The growing size of a positive potential cell is 
equivalent to the motions of the shocks bounding the cell. This can be seen by comparing the innermost shock 
pairs on the right-hand side of the center of the simulation region in Figures 4k and 4 i . The motion of shocks 
enables the positive potential cells to merge together. This elevates the potential of the outer region of the flux tube 
with respect to its ends, as well as with respect to the central region where the potential is negative in a relatively 
extended region. 

We see from Figures 4d to 4f and 4j to 4^ that as the innermost shocks bordering the central 
counterstreaming move inward, the size of the central counterstreaming region decreases. However, this does not 
continue indefinitely. The evolutions of the ion and electron phase space and that of the potential distribution, for 
times later than that shown in Figure 4, are given in Figure 9. Figures 9a to 9c for the ion phase space and the 
corresponding potential distribution in Figures 9g to 9i show that for i > 2200 , the central region of ion 
counterstreaming expands again as the innermost shocks move outward. The corresponding electron phase space 
plots in Figures 9m to 9r show that in the central negative potential region the electrons' velocity spread shrinks. 
Why do the innermost shocks move outward again? The answer to this question lies in the change in the plasma 
influx into the tube caused by the elevation of the positive potential in its outer regions due to the merging of the 
positive potential cells as shown in Figures 4k and 4 1 . The elevated potential chokes the plasma flow into the flux 
tube. When the ion flux is reduced, the inner shocks begin to move outward to maintain the continuity of the ion 
flux through the shocks. This phenomenon of shock motion is quite similar to the motion of double layers seen in 
numerical simulations [Singh and Schunk , 1982] and laboratory experiments [Iizuka et al, 1983]; in these studies 
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it was found that when ion flux to an existing double layer is reduced, the double layer remains stable in its shape, 
but it moves toward the ion source to keep the ion flux through the double layer at a constant value. The shock and 
the double layer are both electrostatic structures and it appears that the shock also remains stable in response to the 
changes in the ion flux by moving in an appropriate direction. For the outward shock motion of the shock in 
Figure 9g to 9j, the shock velocity is found to be from t — 2200 to 2800. At later times the i-i 

instability again inflicts the counterstreaming ion beam as seen from the plots for / > 2800 in Figure 9. This 
mixes the ion beams in the outer region of the counterstreaming, but leaving the ions counterstreaming in the 
central region. It appears that the above process of i-i instability in the outer region, formation of positive potential 
cells and their merging, shock formations and their motion may repeat a few times before the flux tube fills with a 
thermal i zed plasma. 

The formation of shocks in the off-equatorial region is predicted by two-stream hydrodynamic models 
[Rasmussen and Schunk, 1988; Singh, 1990], But these shocks form due to the slowing down of the ion beams by 
polarization electric fields and/or Coulomb collisions. A semikinetic model of interhemispheric plasma flow does 
not show shock formation [Wilson et a!., 1992], Recently Singh et al, [1994] found that in a semikinetic model, 
which treats electrons by assuming a Boltzmann distribution, the condition on the electron to ion temperature for 
the ion-ion instability is not met. Our small-scale simulation discussed here show that temperature conditions can 
be satisfied by the electron heating due to the electron-ion instability which is driven by relatively fast ion beams 
and occurs during the early stage of the refilling. When electron dynamics is simplified by assuming the 
Boltz man n law, this instability is lost along with the associated consequences of electron heating, ion-ion 
instability and shock formation. 

It is important to point out that dynamics of the shocks seen in two-stream hydrodynamic models is quite 
diff erent from that seen in Run-2. In the former case a shock forms in each stream in the off-equatorial region and 
subsequently it moves upward to the equator and then downward to the ionospheric boundary. Such shock motions 
affect the transition from supersonic to subsonic flow for each plasma stream [Singh, 1991a], The resulting 
subsonic flows continue to counterstream indefinitely. On the other hand, Run-2 shows that the counterstreaming 
in the off-equatorial region is rather quickly thermalized by the ion-ion instability and it persists over a longer time 
only in the central region. However, as seen from Figures 4 and 9 the spatial extent of the central 
counterstreaming is quite dynamic. Furthermore, we find that the i-i instability should lead to the formation of 
several shock pairs evolving from the positive potential cells created by the instability. The evolution and merging 
of such shock pairs has a profound effect on the state of the plasma in the central (equatorial) region. 

Flux Tube Filling: The ion-ion instability and its nonlinear evolution affect the filling of the flux tube in several 
ways, which we discuss here. Figures 10a and 10b show the comparison of the ion density distributions from Run- 
1 and Run-2 at 7 = 2000 and 3600, respectively. Figure 10a shows that compared to Run-1, in Run-2 the filling 
is depressed in the central region while it is enhanced in the off-central region at / = 2000 . The explanation for 
this feature of the filling lies in the potential distribution as shown in Figures 4j to At . The elevated potentials in 
the off-central region inhibit the flow of plasma into the central region; this keeps the ion density depressed there 
and therefore the plasma entering the flux tube accumulates in the outer regions of the flux tube. We find that 
until about t = 2000 the net plasma contents of the flux tube in Run-1 and Run-2 are approximately the same. 
This is seen from Figure 1 1, in which the total number of computer ions accumulated in the flux tube is plotted as 
a function of time for the two runs. 

Figure 10b shows that eventually the ion densities in the central region of the flux tube in the two runs 
become nearly equal. This is due to the fact that the potential barriers for the ions in the outer region of the flux 
tube reduces as seen by comparing the potential distributions in Figures 4j, 4k and 4/ for t < 2000 with those 
in Figures 9g to 9i for t > 2000. However, the densities in the off-central region in Run-2 remains generally 
larger than that in Run-1. This yields a larger plasma content in the former run than that in the latter. Figure 11 
clearly shows this for t > 2000. The enhanced refilling in Run-2 is attributed to the trapping of ions and 
electrons in the vortices set up by the ion-ion instability. 
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4. Conclusions and Discussion 

The main conclusions of this paper are the following: 

1. As the magnetic flux tube fills with a plasma expanding from plasma sources located at its end, the ion-beams 
driven instabilities significantly affect (a) the distribution and dynamics of the plasma in the tube, and (b) the 
accumulation of plasma in it. 

2. Two types of instabilities are found to occur: the instability caused by electron-ion interaction is kinetic in 
nature and occurs when the ion beams are relatively fast, and electron dynamics is not simplified by the 
assumption of the Boltzmann law. The ion-ion instability occurs for comparatively slower ion beams. 

3. In the simulation with the Boltzmann electron, the electron-ion instability is suppressed and the ion beams 
remain too fast with respect to the ion-acoustic speed to excite ion-ion waves. However, as the beams slow 
down by the filling of the flux tube, they tend toward the marginal instability, especially in the outer region of 
the flux tube, where the plasma fluctuation is enhanced. This causes an anomalous scattering of the beam 
ions in the outer region. 

4. In simulation with the kinetic treatment of the electrons, e-i instability occurs when the beams are relatively 
fast and it heats the electrons. In addition, some electron energization occurs when they are accelerated by the 
inward pointing electric fields set up by the counterstreaming plasma expansion. When electrons are 
sufficiently heated, increasing the ion-acoustic speed so that V^ r <. 2 C s , the ion-ion instability occurs. 

5. The ion-ion instability is effective in coupling the counterstreaming ion beams in the region where instability 
conditions are met. This sets up potential structures including electrostatic shocks. The central 
counterstreaming is seen to spatially shrink and expand as the ion-ion instability and its nonlinear features 
affecting the plasma evolve. 

6. The trapping of plasma in the potential structures enhances the flux tube refilling and affects the distribution 
of the plasma in it. 

7. In the central region of the flax tube where the magnetic field is minimum, the fast ion beams continue to 
counterstream over a relatively long time. 

Following the work of Banks et al [1971], now there are several studies on shock formation including 
large-scale hydrodynamic [Singh et al , 1986; Gutter and Gombosi, 1991; Rasmussen and Schunk , 1988] and 
kinetic models [Wilson et al , 19992; Singh et al, 1994] and small-scale simulations [Singh et al , 1986; Singh, 
1988]. The usual picture emerging from such studies is the formation of a pair of shocks in the equatorial region 
or in the off-equatorial region, one in each hemisphere. The simulation presented here suggests that when i-i 
instability occurs, it may create numerous shock pairs in the off-equatorial region. Their merging and dynamics 
may introduce both spatial structures in the plasma distribution and temporal features in the evolution of the 
plasma state, both in the equatorial and off-equatorial regions. 

In view of the above results, a pertinent question is whether they are relevant to the refilling of 
plasmaspheric flux tube. Since the ion beams are inherent to the early stage refilling, the answer to this question is 
solidly yes. In view of the crudeness in the temporal and spatial resolution of the large-scale refilling models, the 
most pertinent question is how to include the effects of the instability in such models. At this juncture the answer 
to this question is not very clear, but we know that it must be done by calculating anomalous plasma transport 
coefficients which can be included in the large-scale treatments. The heating of electrons and ions by the 
instability can be included through an anomalous collision frequency like that for ion-cyclotron instability in the 
auroral plasma [Ganguli and Palmadesso , 1987], but how to include the effect of ion and electron trapping in 
potential structures set up by the ion-ion instability remains a challenge. 

Acknowledgment: This work was supported by NASA grant NAGW-2128 made to the University of Alabama in 
Huntsville. 


8 


References 


Banks, P. M„ A. F. Nagy, and W. I. Axford, Dynamical behavior of thermal protons in the mid-latitude ionosphere 
and magnetosphere, Planet . Space Sci 19, 1053, 1971. 

Baker, D. R., Nonlinear development of the two ion beam instability, Phys. Fluids, 1730, 1973. 

Forslund, D. W., and C. R. Shonk, Formation of electrostatic collisionless shocks, Phys. Rev. Lett., 25, 1699, 1970. 

Fried, B. D., and A. Y. Wong, Stability limits for longitudinal waves in ion-beam-plasma interactions, Phys. 
Fluids, 9, 1084, 1966. 

Fujita, T., T. Ohnuma and S. Adachi, Self-oscillations excited by two stream ion-ion instability. Plasma Phys., 19, 
875, 1977. 

Ganguli, S. B„ and P. J. Palmadesso, Plasma transport in the auroral return current region, J. Geophys. Res., 92, 
8673,’ 1987. 

Gresillon, D„ and F. Doveil, Normal modes in the ion-beam-plasma system, Phys. Rev. Lett., 34, 77, 1975. 

Guiter, S. M„ and T. I. Gombosi, The role of high-speed plasma flows in plasmasphenc refilling, J. Geophys. Res., 
95, 10,427, 1990. 

Hasegawa, A., Plasma Instability and Nonlinear Effects, Springer-Verlag, Berlin, Chapter 4, p. 160, 1975. 

Ichimaru, S„ Basics of Plasma Physics, The Benjamin-Cummings Publishing Company, Inc, Reading, MA, 
Chapt. 11,’ 1973. 

Iizuka, S., P. Michelsen, J. J. Rasmussen, R. Schritwieser, R. Hatakeyama, K. Saeki, and N. Sato, Double layer 
dynamics in a collisionless magnetoplasma, RJSO -M-2414, 1983. 

Khazanov, G. V., M. A. Kuen, Yu V. Konikov, and I. M. Sidorov Simulation of ionosphere-plasmasphere 
coupling taking into account ion inertia and temperature anisotropy, Planet. Space Sc,., 32, 585, 1984 

Lin, J., J. L. Horwitz, G. R. Wilson, C. W. Ho, and D. G. Brown, A semikinetic model for early stage refilling: 2. 
Effects of wave-particle interactions, J. Geophys. Res., 97, 1121, 1992. 

McKee, C. F„ Simulation of counterstreaming plasmas with application to collisionless electrostatic shocks, Phys. 
Rev.' Lett., ’l8, 990, 1970. 

Rasmussen, C. E„ and R. W. Schunk, Multistream hydrodynamic modeling of interhemisphenc plasma flow, J. 
Geophys. Res., 93, 14,557, 1988. 

Singh, N„ Comment on "Multistream hydrodynamic modeling of interhemisphenc plasma flow, by C. E. 
Rasmussen and R. W. Schunk, J. Geophys. Res., 95, 17, 272, 1990. 

Singh, N ., Ion-Electron instability of ion-beam-plasma systems, Physics Letters, 67A, p. 372, 1978. 

Singh, N., Role of ion temperature anisotropy in multistage refilling of the outer plasmasphere, Geophys. Res. 
Lett., 18, 817-820, 1991a. 

Singh, N„ Role of electromagnetic noise in the interhemisphenc plasma exchange along closed field lines, 
Advances in Space Res., 11, 9 pp. (9)51 -(9)54, 1991b. 


9 


Singh, N., Refilling of a plasmaspheric flux tube: Microscopic plasma processes, in modeling magnetospheric 
plasma, ed. T. E. Moore and H. H. Waite, Geophys. Monograph , 44, Amer. Geophys. Union, Washington, DC., 
87, 1988.. 

Singh, N., and C. B. Chan, Effects of equatorially trapped ions on refilling of the plasmasphere, J . Geophys. Res., 
97, 1167, 1992. 

Singh, N., and K.S. Hwang, Perpendicular ion heating effects on the refilling of the outer plasmasphere, J. 
Geophys. Res., 92, 13513, 1987. 

Singh, N., W. J. Raitt and F. Yasuhara, Low energy ion distribution functions on a magnetically quiet day at 
geostationary altitude (L=7), J. Geophys. Res., 87, 681, 1982. 

Singh, N., and R. W. Schunk, Dynamical features of moving double layers, J. Geophys. Res., 87, 3561, 1982. 

Singh, N., and R. W. Schunk, Numerical simulations of counterstreaming plasmas and their relevance to 
interhemispheric flow, J. Geophys. Res., 88, 7867, 1983. 

Singh, N., R. W. Schunk, and H. Thiemann, Temporal features of the refilling of a plasmaspheric flux tube, J. 
Geophys. Res., 91 , 13433, 1986. 

Singh, N., and D. G. Torr, Effects of ion temperature anisotropy on the interhemispheric plasma transport during 
plasmaspheric refilling, Geophys. Res. Lett., 17, 925, 1990. 

Singh, N., G. R. Wilson, J. L. Horwitz, Comparison of hydrodynamic and semikinetic treatments for a plasma flow 
along closed field lines, J. Geophys. Res., in press, 1994. 

Wilson, G. R., J. L. Horwitz, and J. Lin, A semikinetic model for early stage plasmasphere refilling, 1, Effects of 
Coulomb collisions, J. Geophys. Res., 97, 1189, 1992. 


10 



FIGURE CAPTIONS 


Figure 1. (a) Geometry of the simulation; the magnetic flux tube and the two plasma reservoirs at the ends of the 
tube are shown. The flux tube fills with plasma when the plasmas from the reservoirs expand into it. (b) Magnetic 
field distribution in the flux tube; the mirror ratio is 10. 

Figure 2. Evolution of ion phase-space in X — V\\ plane from Run-1. 

Figure 3. Plasma density distributions at some selected times from Run-1. 

Figure 4. (a) to (f) Evolution of ion phase-space in the X - V\\ plane from Run-2. Note the coupling between 
the ion beams due to ion-ion instability for / > 1400. (g) to (?) Potential distribution corresponding to the 
phase-space evolution in Figures 4a to 4f. The arrows in panels j, k and ? show the evolution of the shock pairs 
from the positive potential cells created by the ion-ion instability, (m) to (r) Evolution of electron phase-space in 
the X-Fu plane. 

Figure 5. Ion parallel velocity distribution functions from Run-1 at some selected times in the region 
4000 < X < 4500, where scattered ions were seen for I > 4000. 

Figure 6. Temporal evolution of density fluctuations at (a) X = 500, (b) X = 1000 and (c) X = 1500. 
Note the increase in the fluctuation amplitude with decreasing X . 

Figure 7. Ion parallel velocity distribution functions at Junes (a) I = 1000, (b) t = 1200and (c) t = 1400. The 
distributions are for the ions in the region 4000 < X <, 4500. The corresponding electron distributions are 
shown in (d) t = 1000, (e) t — 1200 and (f) / = 1400. 

Figure 8. Temporal evolution of (a) Electron and ion densities, (b) Electric potential <f>, all at X = 2500. (c) 
Frequency spectrum of the potential shown in (b). The dotted curve in (b) shows the average potential when the 
fluctuations are averaged out. 

Figure 9. Same as Figure 4, but the plots are for t > 2000. Note the temporal variation in the spatial extent of 
the central counterstreaming and the associated charges in the potential distribution. 

Figure 10. Comparison of the density distributions in the flux tube from Run-1 and Run-2 at (a) t = 2000 and 
(b) / = 3600. 

Figure 11. Comparison of the flux tube contents from Run-1 and Run-2. Note the enhanced filling in Run-2. 
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Abstract: 

Equatorially trapped hot plasmas are a common feature of the outer plasmasphere, where flux tube 
refilling with cold ionospheric plasma occurs after magnetic storms. The role of the hot plasma consisting of hot 
anisotropic ions and isotropic warm electrons in the refilling process is examined by means of numerical 
simulations using a one-dimensional particle-in-cell code. Simulations are performed on the filling of an artificial 
flux tube having a minimum magnetic field at its center. We have performed two types of simulations; in one type, 
called here Run-A, we allowed cold plasmas to flow into a centrally trapped hot plasma consisting of warm 
isotropic electrons and hot anisotropic ions with perpendicular temperature 7^ > 7j], the parallel temperature. 
Run-A reveals a variety of plasma processes relevant to the plasmaspheric refilling affected by the presence of a hot 
plasma, including formation of propagating electrostatic shocks, intrinsically unstable plasma distribution 
functions produced by the mixing of hot and cold plasmas, weak downward electric fields supported by an 
extended potential distribution in the relatively late stage of the evolution of plasma in the flux tube, and an 
enhanced flux tube filling. In the other type of simulation first a cold plasma flow was allowed to set up in the 
flow, then a hot plasma consisting of the isotropic electrons and anisotropic ions (7^ > 7\\) was suddenly injected 
into the central region of the flux tube. In this case the main distinguishing feature was the formation of relatively 
stable shocks near the mirror points of the centrally trapped hot plasma. The shocks were found to be standing, 
unlike in the previous type of simulation. Since the standing shocks form near the effective mirror points of the 
centrally trapped hot ions, they are called mirror shocks to contrast them from the moving electrostatic shocks seen 
in Run-A. The stability of the standing shocks was found to increase with the decreasing temperature of the warm 
electrons injected with the hot plasma. Wherever possible, similarities between the results from the simulations 
and those from observational data are pointed. 

1. Introduction 

In a companion paper [Singh and Leung, 1994], hereafter referred to as Paper 1, we have studied the role 
of ion-beam driven plasma instability in filling a magnetic flux tube with a cold plasma. In this paper, we extend 
that study by including a hot plasma population trapped in the central region of the flux tube where the magnetic 
field is minimum as shown in Figures la and lb. This study is motivated by the fact that the outer region of the 
plasmasphere is always populated by a hot plasma of plasmasheet and the ring current origins [/toax et al, 1982; 
Deforest et al , 1971]. Furthermore, at times a warm plasma population is created when cold ions in the equatorial 
region are transversely heated by waves generated by the hot plasma population [Olsen et al , 1987]. How do these 
warm and/or hot plasma populations affect the flow of cold plasma and the refilling of the plasmasphere? The 
purpose of this simulation study is to develop an understanding of the plasma processes driven by hot-cold plasma 
interactions, which are relevant to the plasmasphere. 

The simulations show that a hot and/or warm plasma consisting of warm electrons and hot anisotropic 
ions with 7^ > 7j] where 7^ and 7j] are the perpendicular and parallel ion temperatures, respectively, have a 
profound effect on the flow of cold plasma and its accumulation in the flux tube. The hot plasma injection creates 
a potential barrier for the ions in the cold plasma flow. When the cold plasma flow comes into contact with the 
equatorially trapped hot plasma and the associated potential barrier, electrostatic shocks form near the contact 
points. The shocks propagate away from the hot plasma thermalizing the plasma behind it and leaving behind an 
extended potential distribution which is effective in trapping the cold ions and in enhancing the filling of the flux 
tube. On the other hand, when the hot plasma is injected into the cold plasma flow with persistent 
counterstreaming ion beams in the central region of the flux tube [Singh and Leung , 1994], standing shocks form 
with a relatively stable potential structure localized in the minimum magnetic field region. The stability of the 
standing shocks increases with decreasing temperatures of the warm isotropic electrons injected with the hot 
anisotropic ions. 

When the hot and cold plasma spatially mix together, unstable plasma distribution functions are created; 
the salient feature of the ion perpendicular velocity distribution is that a ring distribution forms around the core of 
cold ions. The minimum energy of the ring distribution approximately corresponds to the potential barrier set up 
by the hot anisotropic ions. Such low energy rings can drive lower-hybrid and/or magnetosonic waves below and 
near the lower-hybrid frequency [Lee and Birdsall , 1979], which can heat the cold ions. This is in contrast to the 
energetic ( > 5 keV) ion rings produced by the injection and transport of ions [itoi/x et al t 1982; Deforest and 
Mcllwain, 1971]; such energetic rings create waves which are too fast to resonate with the cold ions and, therefore, 
are ineffective in heating them. 
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The simulation model used in this paper is the same as discussed in section 2 of Paper 1. However, in the 
simulations described here we include a hot plasma in the central region 2300 < X < 2700 (Figure la). The 
properties of the hot plasma and the timing of their injection are described later at appropriate places. Some initial 
results from the present study were reported earlier [Singh, 1993], 

The rest of the paper is organized as follows. In section 2 we describe the results from the simulation in 
which hot plasma is injected initially at time t = 0. In section 3 we describe the results from the simulation in 
which hot plasma is injected at a later time when the cold plasma flow is already set up in the flux tube. The 
properties of the velocity distribution functions of electrons and ions as seen from the runs are described in section 
4. Section 5 deals with the filling of the flax tube with cold plasma. The conclusions of the paper are given in 
section 6. 

2. Cold Plasma Flowing into a Hot Plasma 

We first discuss the results from a simulation in which cold plasmas expanding from the plasma 
reservoirs in Figure la come into contact with the hot plasma trapped in the central region of the flux tube where 
the magnetic field minimizes (Figure lb). Hereafter we refer to this simulation as Run-A. The situation 
considered in this simulation corresponds to that of the fast ion streams flowing into the equatorially trapped hot 
plasma during the early stage of the plasmaspheric refilling. Some preliminary results from this run were 
presented in a previous paper [Singh, 1993], The hot plasma has the following properties; 

T/l = 2 Tff - 1 800 T 0 and Tjj_ = Tj^ = T e = 1 0 T 0 , where T a is the cold plasma temperature in the plasma 

reservoirs. When the hot plasma is the only plasma in the flux tube, the difference between the temperature 
anisotropy of electrons and ions sets up electric fields pointing away from the central region and the corresponding 
potential distribution is given by [Whipple, 1977] 

■Kx)^(kT»U)[\+T«IT,r'tnr (1) 

where r= (Tff / TJj^) (1 - B m l B(x)) + B m / B(x), B m is the minimum magnetic field, and B(x) is the 
magnetic field at the point where the potential is <f>(x) . 

When B m / B « 1 and T e « Tfl , equation ( 1) is simplified to <f>(x ) = (kT e / e ) in ( / Tff) . For 

T e = \0T o , and T^/Tff =2, the potential difference between the point of minimum magnetic field and the 
point where the magnetic field is B(x) is about lkT 0 te. This maximum potential drop is expected without 
any cold plasma in the flux tube. 

Large-scale refilling models [Singh et al, 1986; Wilson et al, 1992] indicate ion flow velocities upto 30 
km/s during the early refilling stage, the corresponding kinetic energy of H + ion beam is < 5 eV, which is 
smaller than the maximum potential energy of the ions estimated above for T 0 - ieV. Thus, the hot plasma is 
expected to have a significant effect on the flow of cold plasma, and hence on the refilling processes. 

We performed several simulations in which hot plasma properties were varied. We found that as long as 
hot ions are anisotropic with 7-^ > Tffi and the warm electrons are isotropic or anisotropic with the reverse 

anisotropy, that is, > Tjj~, the basic effects of the hot plasma on the cold plasma flows are the same. 

Therefore, we describe here only the results from the simulation ran with the hot plasma properties as mentioned 
above. 

Figures 2a to 2i shows the evolution of the ion phase-space in the X -V~\\ plane. Figure 2a for 
t - 200 shows the expanding ion beams from the plasma reservoirs at X = 0 and X - d . For showing some 
details of the cold ion phase space, the vertical scale in Figure 2 is expanded and it does not fully show the velocity 
space of the hot ions. At 7 = 400 (Figure 2b), the ion beams begin to come into contact with the hot plasma. 
The plot at l = 600 (Figure 2c) shows that the incoming ion beams are being reflected. This reflection is 
expected because the potentials set up by the hot plasma provide sufficiently large barriers for the approaching ion 
beams. The evolution of the potential distribution corresponding to that of the phase space (in Figure 2) is shown 
in Figure 3a to 3p. An important feature of the potential distribution is the occurrence of a potential pulse 
coinciding with the central hot plasma for 7 < 1 000 . The peak potential at 7 = 200 is 1 5 kT a / e , which is 
nearly double the potential predicted by equation (1). The development of the large potential is attributed to the 
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space charges developing in the system. In contrast, equation (1) is based on the assumption of a quasineutral 
plasma. The expansion of the cold plasma compresses the potential pulse as seen from the plots in Figure 3 for 
200 <t < 800, after which the edges of the pulse steepen into shock like structures and the pulse broadens again 
and slowly diminishes in amplitude reaching an asymptotic value of about SkT a /e at t = 1400. The 
broadening of the potential pulse is associated with the outward propagation of the shocks. The evolution of the 
potential distribution at later times is affected by the shock propagation and, outside the central region between the 
two sWWc it is also affected by the ion-ion instabilities as discussed in Paper 1. The ion-ion instability creates 
vortices in the X - V A \ phase space as seen at t = 1 600 in Figure 2g and the corresponding perturbations in the 
potential distribution in Figure 3g. 

It is useful to discuss the mechanism of the shock formation. When counterstreaming ion beams begin to 
interpenetrate, a pair of electrostatic shocks form due to the coupling between the ion beams affected by the ion-ion 
instability [Forslund and Shonk, 1970; Singh, 1990], As we discussed in Paper 1, in the absence of the central hot 
plasma, the ion beams continue to counterstream indefinitely because their relative beam velocity V rb is too fast 

to satisfy the conditions for the ion-ion instability which are V rb <2 C 3 and T e >3Tj where C s is the 
ion-acoustic speed and T e and T t are the electron and ion temperatures; the instability in the central region 
does not occur even though the latter condition on temperature is satisfied due to the electron heating by the 
electron-ion instability. In the present situation with the central hot plasma, the counterstreaming on both sides of 
the hot plasma is created by the reflection of the ion beams by the central potential barrier. The relative ion beam 
velocities between the counterstreaming ion beams on both sides of the potential pulse near the reflection point are 
0.5, 0.4, 0.32, and 0.26 at I = 600, 800, 1000 and 1200, respectively. We recall that central hot plasma 
consists of a warm electron population with temperature T e = 107 ^, which yields an ion-acoustic speed 
U = 0. 16. Thus the condition for the ion-ion instability near the reflection point is satisfied for I > 1000. It is 
important to point out that the shocks form just outside the region where hot ions are trapped; outside this region 
the cold ions dominate and satisfy the temperature condition mentioned above. The instability couples and mixes 
the counterstreaming ion beams forming a "mixed" plasma "tab" developing near the shocks, one on each side of 
the central hot plasma (Figures 3e and 3f); the plasma tabs appear in the phase space plots as the darkest areas 
between the hot plasma and the counterstreaming ion beams. The extension of the tabs of mixed plasma along 
with outward propagation of the shocks is clearly seen by comparing the phase space plots in Figure 2 for 
t - 1200 to 2400. The outward propagation of the electrostatic shocks are indicated by the arrows in Figure 
2f to 2i and Figures 3f to it . After their formation at about 7 = 1200, shocks move away from the central 
region with a velocity V sh =0.05V le = F„ , where V te and V ti are the electron and ion thermal speeds in 

the cold plasma, respectively. 

As the shocks propagate out, ion-ion instability occurs in the regions ahead of the shocks as seen at 
7 = 1600, 2000 and 2400 in Figure 2, the instability creates vortexes in the phase space. The combined effects 
of the shock propagation and the ion-ion instability mix the counterstreaming ion beams nearly all along the flux 
tube as seen from Figure 2 for 7 > 2800. This mixing of counterstreaming ion beams all along the flux tube is 
quite different from what we saw in Paper 1 that the counterstreaming persists in the central region (see Figures 4a 
to4f and 9a to 9f in Paper 1). 

The evolution of the electron phase space, corresponding to those of the ion phase space in Figure 2, and 
the potential distributions in Figure 3, is shown in Figure 4. The prominent features of the electrons before 
t < 400 are the two spatially-separated electron populations; the centrally trapped electrons in the potential well 
shown in Figure 3 and the electrons expanding into the flux tube from the plasma reservoirs at X - 0 and 
X - d . For 7 > 1200, we see the signature of the shock pair in the electron phase space; the centrally trapped 
electrons develop sharp edges like the potential distribution (Figure 3) and the ion phase space (Figure 2), all 
marked by arrows. Along with the outward propagating shocks, the centrally trapped electrons expand sway from 
the central region. Even after the shocks have propagated out to the boundaries, the electron phase space shows a 
gradual "bulging" from the ends to the central region, in agreement with the potential distribution in Figure 3 for 
t > 2400 when the potential in the central region remains positive and the potential distribution has gentle 
gradients giving electric fields pointing towards the ends. 

Figures 2 and 4 for 7 < 800 show that as the cold plasmas expand to the central hot plasma, ion 

beams are reflected while the electron gas is accelerated toward the center. Therefore, in the early stage the cold 
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ion beams are initially excluded from the central region, but eventually the electrons pull them into the hot plasma 
region. When the cold electrons are accelerated by the central potential pulse, the additional negative charge 
lowers the peak potential after t > 600 as seen from Figure 3. This allows the cold ions to penetrate into the hot 
plasma region. The penetration is clearly seen from the ion phase plots in the x-v ± plane as shown in Figure 
5. As seen from Figure 2, at times prior to ( = 400 the hot ions trapped in the central region are separated from 
the cold ions in the configuration space. Figure 5a shows the situation when the cold ion beams come into contact 
with the centrally trapped hot plasma, having perpendicular velocities upto about Fj_ = 9 . As the cold plasmas 
expand and the ions reflect from the potential barrier, there are accumulations of plasmas near the reflection points 
forming fronts which move towards the midpoint of the simulation region as seen from the plots for 
600 < / < 1 600 in Figure 5. By the time t = 1 600, the cold ions have punched through the potential barrier, 
and the further evolution of the X - V x phase space does not reveal any new feature. The X - V ± phase space 
plot at the end of the simulation run is shown in Figure 5d. Later we discuss the perpendicular velocity 
distribution function of the ions consisting of the cold and hot populations as seen from this figure. The evolution 
of the electron phase space in X — Vj_ does not reveal any noteworthy feature, except that, like the ions, electrons 
have a warm population trapped in the central region along with the cold electrons. 

Evolution of the plasma density distribution in the flux tube is shown in Figure 6. The density 
distribution at / = 200 shows that the cold plasma is well separated from the centrally trapped hot plasma. The 
maximum normalized density of the hot plasma is about 0.05. As the cold plasmas expand, eventually the density 
maximum at thejnidplane (X = 2500) disappears and a density minimum occurs there. In our simulation, this 
occurs at about t = 800. In a recent paper Olsen et al [1994] suggested that the equatorial density distribution of 
the trapped ions with 7]_ > 7[j depends on the electron temperature anisotropy, specifically they suggested that for 
isotropic electrons the trapped ion density distribution has a density maximum at the equator. On the other hand, 
when T e \\ > T el the density distribution has minimum at the equator. Our simulation shows that the occurrence 

of the density maximum or minimum depends on the stage of the mixing of the cold and the hot plasmas. When 
the hot plasma density dominates, the trapped ion density profile has a density maximum where the magnetic field 
minimizes. When the cold plasma flow begins to dominate the hot plasma density, the maximum disappears and a 
density minimum appears in its place. 

The shocks seen in the ion phase space (Figure 2) and the potential distribution (Figure 3) are also seen 
from the density distributions for t = 2000 in Figure 6; the shocks are indicated by the arrows. Figure 6 shows 
the filling of the flux tube with the cold plasma flowing into it from the plasma reservoirs. Later we discuss the 
filling in more detail. 

3. Hot Plasma Injected into Cold Plasma Flow 

In this section we describe the modification in the cold plasma flow already set up in the flux tube, when a 
hot plasma, like the one described in the previous section, is suddenly injected into the central region of the 
simulation system. Hereafter we refer to this simulation as Run-B. The situation considered in this simulation 
corresponds to the injection of the hot plasma in a relatively late stage of the plasmaspheric refilling. We have 
already discussed in Paper 1 the evolution of cold plasma flow without a hot plasma. So far as the ion phase space 
is concerned, we saw in Paper 1 that after / > 2000 the main feature of the flow is that the ion instability mixes 
the ion beams in the outer region of the flux tube while they continue to counterstream in the central region. We 
also noted that the potential distribution in the central region is negative. The state of the cold plasma just before 
the injection of the hot plasma (/ = 2000) is shown in Figure 7; Figures 7a and 7b show the ion phase space in 
X - V A \ and X - planes, respectively; the corresponding electron phase space plots are shown in Figures 7c 
and 7d; and the potential distribution in Figure 7e. The injection of the hot plasma modifies the cold plasma flow 
in several ways, which are discussed below. 

The evolution of the ion phase space in J- V,-\\ plane after the injection is shown in Figure 8. The 
obvious features of the phase space, as seen from Figure 8, are (1) the coupling between the counterstreaming ion 
beams in the central region where the hot plasma is trapped, and (2) generation of counterstreaming of ion beams 
progressively increasing in extent on both sides of the hot plasma; the extension of the counterstreaming regions is 
marked by arrows in Figures 8a to 8d. However, this newly generated counterstreaming does not last too long; the 
ion-ion instability couples the ion beams and mixes them as seen from the plots in Figures 8d, 8e and 8f. 
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The coupling of the ion beams in the region of the hot plasma is caused by the electric fields generated by 
the hot plasma injection. Figure 9 shows the evolution of the potential distribution after the injection. Note that 
just before the injection the potential distribution has a broad minimum in the central region (Figure 8e). The plot 
at t = 2200 in Figure 9a shows that the hot plasma creates a potential pulse at the center. With increasing 
time, the potential pulse broadens and its edges steepen somewhat. Furthermore, the negative potential regions on 
both’ sides of the potential pulse grow in size from t = 2200 to / = 2800 as indicated by the arrows. The 
regions with negative potentials contain the counterstreaming ion beams shown in Figure 8 for the same time 
interval. We point out here that the growing size of the negative potential region along with the ion 
counterstreaming in it was seen in Run-2 of Paper 1 without the hot plasma as shown in Figure 9 of that paper. 
Therefore, the extension in the size of the counterstreaming is not caused by the not plasma injection, but it is 
caused by the outward motion of the shock-like structures marked by the arrows in Figure 7e due to the decrease of 
the ion flux into the central region by the vortices set up in the outer regions of the flux tube (Figures 7a and 7e). 
This point is discussed in greater detail in Paper 1 . The fluctuation in the potential distribution on both sides of the 
pulse (Figures 9e and 9f) and the associated vortices in the phase space (Figures 8e and 8f) are caused by the ion- 
ion instability, which is eventually seen to affect the potential pulse at late times t > 3400, but the core of the 
pulse is found to remain stable with an amplitude of about 2 . 5 kT 0 / e . 

Injection of Hot Plasma with Colder Electrons: We repeated the above simulation (Run-B) by injecting a hot 
plasma with electrons colder than that in Run-B; we chose T e -3.3T 0 in contrast to 10 7^ in Run-B. We refer 
to this new run as Run-C. As before the hot plasma was injected at t — 2000 and the basic features of the 
evolution of the plasma and potential prior to t = 3000 were nearly the same in Run-B and Run-C. However, in 
Run-C with the colder electrons the central potential pulse remains quite stable without being severely affected by 
the ion-ion instability. This is shown in Figure 10; Figures 10a, 10b and 10c show the potential distributions at 
t - 3200, 3600 and 4000; Figures 10d, lOe and lOf show the corresponding ion phase space in X - plane; 
and similarly Figures 10g, lOh and lOi show the electron phase space. Note that Figure 10 shows the pulse 
evolution over a longer time than that shown in Figure 9 for Run-B with the warmer electrons having T e = 1 0 T 0 . 
The potential distributions show that the stable potential pulse is flanked by negative potential regions in which 
counterstreaming ions persist. In the central positive peak region electrons are trapped electrostatically (Figures 
10g, lOh and lOi), while ions are magnetically trapped in the same region. In the neighboring negative potential 
regions ions are trapped. In view of the alternate electrostatic trapping of the electrons and ions in their respective 
potential wells, the central peak is an "electron hole" and the negative potential wells on both sides of the potential 
pulse are "ion holes." This combination of the central electron and ion holes is remarkably stable at least over a 
time period At - 20,000 which is about 160 ion plasma periods. For a local plasma density of 10 cm -3 , we 
expect that the pulse should be stable for a time period of at least a few hundreds of milliseconds, and it should be 
easily observable from satellites. Olsen et al, [1987] have reported examples of equatorial plasma showing a sharp 
transition between counterstreaming field-aligned flows and an equatorially trapped highly anisotropic warm ion 
population. The central distribution of plasma and the potential, as shown in Figure 10, are qualitatively similar to 
the situation reported by Olsen et al [1987], although at a quite different spatial scale. Our simulations show that 
the transition between the field-aligned counterstreaming and the trapped ions occur near the effective minor 
points of the latter ions. Therefore, we call the transitions mirror shocks to contrast them from the propagating 
electrostatic shocks seen in Run-A. 

4. Plasma Distribution Functions 

We discuss here some interesting features of the electron and ion velocity distribution functions produced 
by the hot plasma injection. We found that the results on the velocity distribution functions from the runs 
discussed here did not differ significantly after the spatial mixing of the hot and cold plasmas. Therefore, in this 
section we primarily discuss the distribution functions as seen from Run-A. 

We examine the distribution functions from two regions of the flux tube; near the ends of the tube, say for 
20 < x < 240, and in the central region covering 2300 < x < 2700. The evolution of the ion distribution 
functions near the ends of the flux tube is shown in Figure 11. At early times (t ~ 200) the distribution shows 
an ion beam flowing into the flux tube. Eventually precipitating ions appear as an ion beam for V\\ < 0; at 
t = 1000 the beam velocity V b » 0. 24 V te = 5F ft . Later on both the inflowing and the precipitating ion beams 
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are seen to broaden progressively with increasing time. Eventually the ion beams for V\\ <0 and fy > 0 begin 
to merge together and by the end of the simulation, the beams have completely merged. The above evolution of the 
ion distribution function is the combined effect of the expansion of the cold plasma into the flux tube, loss of some 
of the hot ions injected into the magnetic trap, reflection of the expanding ion beams from the electrostatic shocks 
(see Figure 2) and the ion-ion instability. A noteworthy feature of Figure 11 is that the precipitating (V\\ < 0) 
ions have energies greater than the inflowing ions. This merely reflects that some of the precipitating ions have 
undergone acceleration by the potential distribution set up by the hot plasma. In the course of the simulation, the 
midpoint (X = 2500) potential varies from about 18 kT a !e to 2.5 kT 0 /e (Figure 2). Thus, the 
precipitating ions cover a spectrum of energy extending to about 15 kT Q / e which corresponds to a maximum 
velocity of about 0. 3 as seen in Figure 1 1 for V\\ < 0 . 

The electron velocity distribution function near the ends of the flux tube is shown in Figure 12. The main 
feature of the electron distribution function is that at early times t = 200, the distribution is symmetric about 
V\\ = 0 and it progressively becomes more and more asymmetric with energized electrons appearing for Vy < 0 . 
The energization appears to occur primarily by the ion-ion instability, which transfers part of the ion beam energy 
to the electrons. The ion-ion instabilities are clearly seen from the vortex formation in the X - phase space 
shown in Figure 2. Electron-ion instability is also likely to heat the electrons, especially during the early stages 
when the fast ion beams begin to precipitate near X = 0 as seen from Figure 1 1 for ( = 10 3 . 

The asymmetry of the ion and electron velocity distribution functions in the relatively late stage of the 
filling of the flux tube, as shown in Figures 11 and 12, suggests that the plasmasphere refilling process in the 
presence of equatorially trapped hot plasma generates a downward heat flow, which may affect the energy budget 
of the topside ionosphere. 

The evolution of the perpendicular velocity distribution functions of the ions and electrons near the ends 
of the flux tube do not shown any interesting features. Examples of such distribution are shown in figures 13a and 
13b, which are approximately the distributions in the cold plasma reservoirs. 

The mixing of the hot and cold plasmas in the central region produces distribution functions having a core 
of cold ions and a halo of hot ions. Such parallel and perpendicular velocity distribution functions are shown in 
Figures 14a and 14b, respectively. (These distribution functions are for the ions in the central region 
(2300 < X < 2700). In the parallel velocity distribution the cold core and the relatively hot ion populations are 
merged, but in the perpendicular velocity distribution, the hot ions are well separated from the cold population and 
they extend from a low velocity of V± = 0. 2 to large velocities corresponding to that of the injected hot plasma. 
The hot ion population in Figure 14b appears as a beam in Vj_. Since in gyrophase these hot ions are randomly 
distributed, the beam is a ring distribution in the perpendicular velocity. The ring distributions of energetic ions 
with energies 5-30 keV have been observed in the equatorial region of the outer plasmasphere and they are known 
to excite fast magnetosonic waves [Perraut et al, 1982], However, the ring distribution seen from the simulation is 
considerably less energetic in the sense that its low energy end begins at an energy determined by the peak 
potential of the positive potential pulse seen in the simulations (Figures 3, 9 and 10). The production of such a 
ring distribution function has a simple explanation. The charged particles trapped inside a magnetic mirror have a 
loss cone distribution extending to zero velocities as long as the pitch angle (a) of the ions at the point of minimum 
magnetic field ( B min ) satisfies the condition a> sin -1 (R~ ]n ) where R = B m&x / 5 min and 5 max is the 
maximum magnetic field in the flux tube (Figure lb). However, the central potential pulse sets up electric fields, 
which accelerate the ions with relatively low perpendicular energies into the loss cone and they precipitate out. 

In view of the amplitude of the potential pulse evolving with time (see Figure 3), the maximum parallel 
energy gain for the ions occurs during the early stage t < 1 000 . For the maximum pulse amplitude of ~<j> 0 = 1 8 
(Figure 3), the parallel velocity of the ions accelerated by the potential pulse is 
^llmin = (2 q (p 0 / ntj = 0. 3V te . If such ions remain trapped inside the magnetic mirror as the potential pulse 
evolves, the minimum perpendicular velocity at the midplane (X = X mid ) is given by 

*imin = (2 q <P 0 / tan a 0 

where 0C o — sin - * (T? - *^ ). In our simulation R = 10 and we estimate that = 0. lV te . Figures 5 and 14 

show that the trapped particles have perpendicular velocities V ± > 0.2V /e . The discrepancies between this value 
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of V ± and the above estimate for V lmin is attributed to the temporal evolution and the fluctuations of the 
potentials in the flux tube. 

An interesting consequence of the low energy ring distribution is that it can excite a broadband waves 
over a frequency range from ion cyclotron frequency to the lower hybrid frequency [Lee and Birdsall, 1979]. 
However, in our one-dimensional simulation the wave generation process is suppressed. If the waves are 
generated, they can heat the cold ions, merging the cold core with the ring distribution. Composite perpendicular 
velocity distribution functions having a cold core and a halo of warm ions have been measured from the Dynamic 
Explorer- 1 Satellite [Olsen et at, 1987], 

The velocity distribution functions of the electrons in the central region are shown in Figures 15a and 15b. 
The parallel velocity distribution (Figure 15a) function shows a flat top distribution, indicating the trapping of the 
electrons in a potential well set up by the positive potential in the central region. The perpendicular distribution in 
Figure 15b shows a cold core along with a ring of warm electrons, which are trapped inside the magnetic mirror by 
virtue of their large perpendicular velocities. The ring in the perpendicular velocity of the electrons have the 
potential of driving electron cyclotron waves [Hasegawa, 1975]. 

5. Flux Tube Refilling 

In Paper 1 we saw that the plasma waves enhance the filling of the flux tube by trapping plasma in 
vortexes set up by the ion-ion instability. We find that the hot plasma injected into the central region of the flux 
tube has a similar effect, although through a different process. We compare the filling of the flux tube as seen from 
the simulations reported in this paper (Run-A and Run-B) and those in Paper 1. Figure 16 shows the comparison 
of the flux tube contents from Run-2 of Paper 1, in which the filling was studied without the hot plasma, and Run- 
A and Run-B of this paper. Curve 1 is the same as in Figure 12 of Paper 1, but it is translated along the vertical 
axis by the initial number of the hot ions injected into the flux tube at 1 = 0 in Run-A. Curves 2 and 3 show the 
evolution of the plasma content from Run-A and Run-B, respectively. The vertical translation of Curves 1 is done 
for the purpose of easy comparison of the cold plasma contents of the flux tube as seen from the different runs. We 
point out that the plasma content is just the total number of ions in the flux tube. Since each particle used in the 
simulation is equivalent to a large number of the real ions, we call the simulation ions macroions or computer ions. 
A comparison of Curve 1 with Curve 2 shows that the initially injected hot plasma facilitates the trapping of more 
cold plasma than that from the run without the hot plasma. In the simulation with the initially injected hot plasma 
(Run-A) enhanced trapping occurs by two processes: (1) the vortices set up by the ion-ion instability like that 
shown by the enhanced filling in Run-2 (Curve 1) over that in Run-1 of Paper 1 (Figure 12 Paper 1); (2) the 
potential pulse and its evolution in an extended potential distribution (Figure 3) provide an additional mechanism 
for the trapping of the cold plasma. When cold plasma flows into the flux tube, the electric fields set up by the 
extended potential distribution slowly retard the flow and increase the ion density in the flux tube. The effect of 
this additional mechanism on the filling is illustrated in Figure 17, in which we compare the density profiles at the 
end of the simulations in Run-2 of Paper 1 and Run-A described in this paper. We notice that the enhanced 
density is confined in the central region from Z = 10 3 to 4X10 3 . One can argue that this enhanced density 
is due to the centrally trapped hot plasma, but the comparison of the total plasma contents from Run-2 and Run-A 
in Figure 16 clearly points out that the cold plasma content is enhanced. Since in Figure 17 the relative density 
enhancement in Run-A occurs over an extended region only inside the central portion of the flux tube, it is obvious 
that an additional amount of cold plasma is also contained in that region. The containment of the additional cold 
plasma is facilitated by the retardation of the inflowing ions by the relatively weak electric fields associated with 
the extended potential distribution in Figure 3 for t > 2400 

In Run-B, in which the hot plasma was injected after the cold plasma flow was well set up in the flux tube, 
an enhanced refilling did not occur. In Figure 16, Curves 2 and 3 show the comparison of the total plasma 
contents from Run-A and Run-B. Curve 3 (solid line) shows the sudden injection of the hot plasma at t = 2000. 
With the injection the total plasma content suddenly jumps from Curve 3 to Curve 1, but the content given by 
Curve 3 (Run-B) remains slightly lower than that given by Curve 2 (Run-A) near the injection time and it 
continues to show a lower value of the plasma content till the end of the simulation. We note that the contents 
shown by Curve 3 and Curve 1 for t > 2000 are quite close and are lower than that given by Curve 2 (Run-A), 
indicating that the late injection of the hot plasma does not contribute to the enhancement in the filling with the 
cold plasma. This is understood by noting that in Run-B the potential peak generated by the hot plasma does not 
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evolve into propagating shocks, creating an extended potential distribution with weak electric fields, which plays a 
key role in retarding the inflowing ions and in enhancing the cold plasma density in Run- A. 

6. Conclusion and Discussion 

The main aim of this paper is to study how a hot plasma trapped inside the minimum magnetic field 
region of a flux tube affects the filling of the tube with a cold plasma. We performed two simulations: (1) by 
injecting hot plasma initially and then allowing the cold plasma to expand into it, and (2) first a cold plasma flow 
was allowed to establish in the flux tube and then the hot plasma was injected. The simulations presented here, 
with the initial and delayed injection of the hot plasma, are illustrative examples of how hot and cold plasmas may 
come into contact. In the case of the outer plasmasphere, the contact may occur in a variety of situations ranging 
from the initial refilling stage with the cold plasma density < 1cm -3 to that of a nearly filled plasmasphere with 
densities of ~ 10 3 c/n“ 3 . The hot plasma in both cases consisted of hot anisotropic ions with 7j_ > Tjj and warm 
isotropic electrons. The simulations reveal new processes driven by the hot-cold plasma interactions. These 
processes are summarized as follows. 

Flow of Cold Plasmas into a Trapped Hot Plasma: 

1) Hot plasma creates a potential barrier for the cold plasma flow. The magnitude of the potential differs from 
that given by theory [ Whipple , 1977] which does not account for the presence of the cold plasma. We show 
that the potential evolves as the hot and cold plasmas interact. 

2) The main feature of the evolution of the potential distribution is the formation and propagation of the 
electrostatic shocks. We saw in Paper 1 that in the absence of the hot plasma, the counterstreaming expansion 
of cold plasma in the flax tube sets up counterstreaming ion beams, which last for a relatively long time in the 
central region. In contrast to this, the potential barrier set up by the hot plasma reflects the incoming plasma 
flows on each side of the centrally trapped hot plasma. In the early stage the cold plasmas do not penetrate 
into the central region of the hot plasma. The reflection of ion beams sets up counterstreaming on each side of 
the hot plasma. The shock formation begins in the reflection region where beams are relatively slow and the 
warm electrons give a relatively large ion-acoustic speed. These two effects create the conditions for the ion- 
ion instability (see Paper 1), and hence the shocks form. 

3) The propagation of the shock extends the positive potentials created by the hot plasma to large distances from 
where the hot plasma is trapped. 

4) As the shocks propagate away from the central region of the flux tube, ion-ion instability occurs in the outer 
regions, as described in Paper 1. The combined effects of the shocks and the vortexes set up by the instability 
create an extended potential distribution, in which central region of the flux tube remains generally positive. 

5) The extended potential distributions help in trapping more cold plasma than that in the run without the hot 
plasma. The trapping occurs by retarding the inflowing ions in the flux tube. Thus, the filling of the flux tube 
with the cold plasma is enhanced by the presence of the equatorially trapped hot plasma. 

6) The plasma distribution functions resulting from the hot-cold plasma interaction show some interesting 
features. We found that the precipitating ion beams and electrons, which are leaving the simulation system, 
are energized by the plasma processes occurring in the flux tube. This suggests that the hot-cold plasma 
interaction my eventually deposit some energy into the topside temperature during the refilling of a 
plasmaspheric flux tube. 

7) The ion velocity distribution function in the central region consists of a cold core and a halo of hot ions. In the 
perpendicular velocity distribution, the hot ions appear as a ring starting at a low perpendicular energy which 
approximately corresponds to the maximum potential barrier set up by the hot plasma population. Such ring 
distributions can generate waves in the frequency band from the ion gyrofrequency to the lower-hybrid 
frequency [Lee and Birdsall 1979], However, such waves cannot be studied by the one-dimensional 
simulations used in the present work. 

8) Even the electron perpendicular velocity distribution shows a ring, which are known to drive electron 
cyclotron waves at frequencies / = (n + f ce , where f ce is the electron cyclotron frequency and n is 
an integer [ Hasegawa , 1975]. 
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Injection of a Hot Plasma into a Cold Plasma Flow: 

When a hot plasma is injected into a cold plasma flow already set up in the flux tube, we found the 
following interesting results: 

1) The hot plasma with anisotropic ions and warm isotropic electrons generates a positive potential pulse, which 
is relatively stable as it remains confined to the central region. The stability of the pulse increases with colder 
electrons in the injected hot plasma. The edges of the pulse appear like a pair of standing shocks with 
counterstreaming ions on both sides of the potential pulse. This situation is qualitatively similar to that found 
in the equ ato rial plasmasphere with sudden transition between counterstreaming ions and the equatorially 
trapped anisotropic ions with 7j_ > 3jJ [Olsen et al, 1987J. 

2) The standing shocks form near the boundaries of the centrally trapped hot ions. Since the ions are trapped by 
the mirror face, the shocks form near the effective mirror points and, therefore, we call them mirror shocks. 
This no m e n cl at ure distinguishes these shocks from the electrostatic shocks created by the electrostatic 
processes of ion-ion instability. 

3) The injection leads to expulsion of cold plasma, reducing the plasma content of the flux tube relative to that in 
the run with the initial injection of hot plasma. 

4) The mixing of hot and cold plasma eventually generates plasma distribution functions similar to those shown 
in Figures 14 and 15 for the run with the initial injection of the hot plasma. 

The above results cannot be directly applied to the problem of plasmaspheric refilling. The purpose of 
this paper is to study the variety of plasma processes driven by hot-cold plasma interactions which occur in the 
outer plasmasphere. Some of the processes, such as the formation of shocks and the creation of the potential 
distributions in the flux tube capable of trapping more plasma are well brought out by the one-dimensional 
gimnlatinn On the other hand, the results presented above indicate the possibility of some processes which need 
multidimensional simulations. For example, the distributions shown in Figure 14 are inherently unstable and they 
can excite waves in the frequency band from proton gyrofrequency to the lower hybrid frequency. Such waves have 
been observed in the equatorial region and have been invoked to explain the heating of the cold ions [Olsen et al, 
1987], However, in the past it has been suspected that the waves are driven by energetic (5-20 keV) proton ring 
distributions. The waves excited by such beams are too fast to resonate with cold ions and therefore they cannot 
heat them. The low energy rings generated by the hot-cold plasma interactions can excite relatively slower waves, 
which could heat the cold ions merging the cold plasma with the hot plasma ring in Figure 14. The merging is 
likely to produce composite distribution functions with two temperatures as reported by Olsen et al [1987], In a 
recent paper Olsen et al [1994] suggested that the density distribution of the trapped ions in the equatorial region 
depends on the electron temperature anisotropy. We found from our simulations that the nature of the density 
distribution in the presence of the trapped ions depends on the evolution of the cold plasma density relative to that 
of the trapped ions. When the latter density dominates a density maximum occurs at the equator, otherwise a 
density minimum occurs when the cold plasma density dominates. This is found to be true when the electron 
distribution has both trapped and field-aligned populations. 
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FIGURE CAPTIONS 


Figure 1. Geometry of the simulation, (a) The flux tube with the cold plasma reservoirs near its ends are shown. 
The hot plasma is injected in the center of the flux tube where the magnetic field is minimum as shown in (b). 

Figure 2. Temporal evolution of ion phase space in the X - plane (Run- A). Note the central hot ions; the 
entire distribution of the hot ions is not shown; it is truncated at Vj\ \ = il so that the cold plasma distribution can 
be adequately resolved. Note the reflection of ions, formation and propagation of electrostatic shocks. 

Figure 3, The temporal evolution of the potential distribution corresponding to that of the ion phase space in 
Figure 2. Note the evolution of the shock from the potential pulse set up by the centrally trapped hot plasma. 

Figure 4. The temporal evolution of the electron phase space in X - V e \\ plane corresponding to that of the ions 
in Figure 2 and the potential distribution in Figure 3. 

Figure 5. Evolution of ion phase space in the X-Fj_ plane (Run-A). Note the presence of the centrally trapped 
hot ions. At early times the expanding cold ions create density fronts when they come into contact with the hot 
plasma as shown by arrows. The fronts eventually punch through the hot plasma mixing the cold and hot plasmas 
spatially. 

Figure 6, Temporal evolution of the total ion density in the flux tube (Run-A). At early times the hot plasma 
creates a density maximum at the center. As the cold plasma fills the tube, the density maximum eventually 
disappears. The signature of the electrostatic shocks in the density distribution is indicated by the arrows on the 
density profile for t = 2000 . 

Figure 7. State of the cold plasma in the flux tube at / = 2000 just before the delayed injection of the hot 
plasma (Run-B). (a) X - ion phase space, (b) X - V iL ion phase space, (c) X - 7~ e \\ electron phase space, 
(d)X-^ x electron phase space, and (e) potential distribution. 

Figure 8. Temporal evolution of the ion phase space in X - f^ji plane (Run-B). Note the central mixing of the 
counterstreaming ion beams shown in Figure 6a in response to the delayed injection of the hot ions. Also 
noteworthy, is the extension of the counterstreaming on both sides of the hot plasma and eventual mixing by the 
ion-ion instability. 

Figure 9. Evolution of the potential distribution in response to the delayed injection of the hot plasma (Run-B). 
Note the appearance of the potential pulse in the central region. 

Figure 10. Stability of the potential pulse and its evolution into standing mirror shocks are shown when the 
temperature of the warm electrons in the hot plasma is reduced from 10 T a in Run-B to 3 T a in Run-C. (a) to (c) 
Potential distributions, (d) to (0 Ion phase space in X - plane, (g) to (i) Electron phase space in X - ^ 
plane from Run-C. 

Figure 11. Evolution of the ion velocity distribution function F i (V\\) near the end of the flux tube ( 
100 < X < 500, Run-A). Note the ions reflected ( < 0) from the central potential pulse created by the hot 

plasma injected at t = 0. Eventually the inflowing and the reflected ions merge together. Asymmetry of the 
distribution function suggests heat flux into the topside ionosphere. 

Figure 12. Same as Figure 1 1, but for the electrons. 
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Figure 13. Perpendicular velocity distributions (Run-A). (a) Ions and (b) Electrons (100 < X < 500) at 
f = 3600. 


Figure 14. Ion velocity distribution functions at / = 3600 in the central region of the flux tube where hot and 
cold ions have mixed, (a) F t ( V \\ ) and (b) (V ±) . Note the presence of the cold core ions and the hot ions 

appearing as a halo in F i ( V \\ ) and a beam in F i (Kj_) . 

Figure 15. Same as Figure 14, but for electrons, (a) F e (V\\) . (b) F e . Note that the electrons trapped with 
the ions have a flat top distribution in F e (V\\) and a beam in F e (V^). 

Figure 16. Comparison of the flax tube plasma contents from Run-2 of Paper 1 and Run-A and Run-B of this 
paper. Note that Curve 1 from Run-2 is vertically translated by adding the number of injected hot ions to its 
plasma content. This easily brings out the additional filling with the cold plasma caused by the processes driven by 
the hot plasma. After the hot plasma is injected at / = 2000 in run-B, the plasma content shown by Curve 3 
closely follows Curve 1, and both remain below Curve 2. 

Figure 17. Comparison of the density distributions from Run-2 in Paper 1 and Run-A with the injection of the hot 
plasmas at t = 0. The densities shown are for t = 3600. 
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